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ABSTRACT 

J 

An  exponential  finite  difference  algorithm,  as  first  presented  by 
Bhattacharya  for  one-dimensional  unsteady  state,  heat  conduction  In 
Cartesian  coordinates,  has  been  extended.  The  finite  difference 
algorithm  developed  was  used  to  solve  the  diffusion  equation  In 
one-dimensional  cylindrical  coordinates  and  applied  to  two-  and 
three-dimensional  problems  In  Cartesian  coordinates.  The  method  was 
also  used  to  solve  nonlinear  partial  differential  equations  In  one 
(Burger's  equation)  and  two  (Boundary  Layer  equations)  dimensional 
Cartesian  coordinates.  Predicted  results  were  compared  to  exact 
solutions  where  available,  or  to  results  obtained  by  other  numerical 
methods.  It  was  found  that  the  exponential  finite  difference  method 
produced  results  that  were  more  accurate  than  those  obtained  by  other 
numerical  methods,  especially  during  the  Initial  transient  portion  of 
the  solution.  Other  applications  made  using  the  exponential  finite 
difference  technique  Included  unsteady  one-dimensional  heat  transfer 
with  temperature  varying  thermal  conductivity  and  the  development  of 
the  temperature  field  In  a  laminar  Couette  flow.  ^ 
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NOMENCLATURE 


aj.b^'C^  Thomas  algorithm  variables 
8  Blot  number 

Cp  material  specific  heat,  J/kg  •  K  (Btu/lbm  •  °F) 

h  convection  heat  transfer  coefficient,  W/M2  •  °C 

(Btu/ft2  •  hr  •  °F) 

1,J,k  nodal  location  In  x,y,  and  z  spatial  coordinate  directions 
respectively 

Jq.J]  Bessel  functions  of  zero  and  first  order 

k  thermal  conductivity,  W/M  •  #C  (Btu/hr  •  #F  •  ft) 

k?  thermal  conductivity  at  1th  position,  nth  time  step, 

W/M  •  °C  (Btu/hr  •  °F  •  ft) 

L  distance  between  plates,  M  (ft) 

M  dimensionless  drive  number 

m  number  of  sub-intervals 

N  number  of  nodes  In  a  spatial  direction 

n  time  step  position  designation 

q  heat  flux,  W/M2  (Btu/hr  •  ft2) 

r  spatial  coordinate;  cylindrical  coordinates,  K  (ft) 

T  temperature,  °C  (°F) 

t  time,  sec 

At  time  between  time  steps  n  and  n  ♦  1 

U  Couette  flow  velocity.  M/s  (ft/s) 

x,y,z  spatial  coordinates,  Cartesian  coordinates,  M  (ft) 

Ax, Ay, A z  distance  between  nodal  positions  In  the  x,y,  and  z  spatial 
directions 

a  thermal  dlffuslvlty,  M2/s  (ft2/s) 

8  rate  of  thermal  conductivity  variation 

y  At/pCp( Ax)2;  (W/M  •  “C)"1  ((Btu/hr  •  °F  •  ft2)-1) 


k  constant 

constant  used  In  exponential  finite  difference  method  with 
temperature  varying  thermal  conductivity 

Ax  finite  difference  operator 

^  rath  eigenvalue  of  Bessel  function 

x?,v?  Thomas  algorithm  variables  dependent  on  time  step  and  spatial 
1  location 

5  amplification  factor 

v  kinematic  viscosity,  m2/s  (ft2/s) 

p  material  density,  kg/M3  (lbm/ft3) 

<p,*,e  separation  variables 

Q  tt  dimensionless  time 

(ax)' 


v  1 


I.  INTRODUCTION 


Partial  differential  equations  have  many  Important  applications  In 
the  fields  of  engineering  and  physics.  Many  exact  solutions  exist 
depending  on,  the  partial  differential  equation,  the  boundary 
conditions,  and  the  number  of  spatial  dimensions  under  consideration 
[1].  Coordinate  systems  other  than  Cartesian,  more  than  one  spatial 
dimension,  and  the  boundary  conditions  all  can  pose  problems  that  are 
either  extremely  difficult  or  Impossible,  to  solve  by  analytical 
methods.  Numerical  methods  thusly  become  the  only  possible  solution 
method  If  the  problem  complexity  Is  not  to  be  compromised.  Typically 
the  ability  of  a  particular  method  to  predict  a  field  variable  Is 
tested  by  numerically  solving  a  problem  for  which  a  known  exact 
solution  Is  available.  The  ability  of  the  method  to  predict  the  exact 
results  Is  a  measure  of  the  confidence  that  can  placed  In  a  solution 
where  no  exact  solution  exists  or  experimental  test  results  are 
unavailable. 

The  objective  of  the  work  to  be  presented  Is  to  extend,  expand, 
and  compare  an  explicit  exponential  finite  difference  technique  first 
proposed  by  Bhattacharya  [2].  To  date  the  method  has  only  been  used 
for  one-dimensional  unsteady-state,  heat  transfer  problems  In 
Cartesian  coordinates. 

The  method  has  been  expanded.  In  this  report,  to  allow  application 
to  a  variety  of  problems.  The  exponential  method  will  be  extended  here 
to  the  case  of  one-dimensional  unsteady  heat  transfer  In  cylindrical 
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coordinates.  Also,  It  was  used  In  two-  and  three-dimensional  unsteady 
heat  condition.  Other  cases  of  Interest  that  were  solved  using  this 
approach  Include  temperature  varying  conductivity  in  one-dimensional 
heat  transfer  and  the  development  of  the  temperature  field  in  laminar 
Couette  flow.  Solutions  of  the  above  cases  were  either  compared  to 
exact  solutions  or  to  results  obtained  by  alternative  numerical 
techniques. 

One  final  application  of  the  exponential  finite  difference 
algorithm  was  made  for  nonlinear  partial  differential  equations. 
Burger's  equation  along  with  the  boundary  layer  equations  are  solved 
using  the  exponential  method.  Thus,  a  demonstration  of  how  to  apply 
the  method  to  nonlinear  problems  Is  described. 

The  results  of  all  the  different  cases  considered  In  this  study 
Indicated  that  the  exponential  technique  produced  results  that  were 
more  accurate  than  those  found  through  other  numerical  techniques.  All 
exponential  computer  codes  and  those  of  the  other  competing  numerical 
analysis  were  run  In  double  precision  on  either  the  IBM-3033  or  the 
Cray  X-MP  mainframes.  The  computer  codes  developed  for  the  exponential 
finite  difference  method  and  other  numerical  techniques  used  for 
comparison  are  contained  In  the  appendix  of  this  report. 


II.  -  ANALYSIS 


The  Exponential  Finite  Difference  Algorithm 
An  explicit  exponential  finite  difference  algorithm  as  first 
derived  by  Bhattacharya  [2]  will  now  be  presented.  The  method  can  be 
applied  to  many  of  the  partial  differential  equations  found  In 
engineering  and  physics.  The  diffusion  equation  as  It  applies  to 
conduction  heat  transfer  will  be  used  In  the  demonstration  that 
follows.  In  reference  [2-3]  the  method  was  derived  for  one-dimensional 
conduction  heat  transfer  In  Cartesian  coordinates.  To  show  how  the 
method  can  be  extended,  a  derivation  parallel  to  the  one  presented  In 
reference  [2]  will  be  made  for  unsteady  state  heat  conduction  In 
two-dimensions.  Equations  of  this  type  are  typically  solved 
numerically  by  a  variety  of  methods  [4]. 

For  two-dimensional  heat  transfer  In  Cartesian  coordinates  with 
constant  material  properties,  the  appropriate  partial  differential 
equation  Is  [5]: 


il 

at 


(i) 


To  Initiate  the  exponential  method  a  product  solution  Is  assumed  and  Is 
written  as: 


T(x,y,t)  =  4>(x)4r(y)e(t)  (2) 


The  Initial  conditions  of  the  problem  are  assumed  to  be 
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Now  taking  the  appropriate  devlatlves  of  Eq.  (2)  with  respect  to  the 
Independent  variables, yields 


dT 


ae 


it  3  *  iit  ; 


32T 


= 


. 

2  ' 


a2T 


ax  ax  ay 

Substituting  Eq.  (4)  Into  Eq.  (1)  produces: 


al* 

ay2 


(4) 


4* 


Dividing  both  sides  of  the  above  by  4*90  gives: 


1  30  (  1 

0  at  B  “|4 


2  2  \ 

LA  .  1  Uf  (  _ 

ax2  *a,2(" 


(5) 


It  can  be  seen  that  the  variables  have  been  separated.  Consequently 
both  sides  of  Eq.  (5)  must  equal  a  constant,  say,  *. 

Now  examining  only  the  left  hand  side  of  Eq.  (5), 


I  ae 
e  at 


-  K 


Multiplying  the  left-hand  side  of  this  equation  by  gives: 

ae 

e**  at  3  "  K 

which  can  be  rewritten  from  Eq.  (4)  as: 

1  aj 

T(x.y.t)  at  x  K 
Direct  Integration  produces: 


T  =  c2  exp  j-  ict j 
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about  a  node  (1,J)  as  [6]: 


Thus  Eq.  (8)  becomes: 
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Now  Eq.  (10)  Is  substituted  to  replace  the  constant  *  In  the 
exponential  of  Eq.  (6).  Making  the  appropriate  substitutions  results 
In: 


T1 J  *  Tin.j  6XP 


aAt 


U  L 


Tn  +  Tn  -  2Tn 

litU _ i-i  J _ 1L1 


(  AX ) ' 


Tn  +  Tn  -  2Tn 

i.J*i  i.J-i _ 1 AA 


(ay)‘ 


If  the  grid  spacing  Is  constant  (Ax  =  Ay),  then  the  above  equation 
can  be  written  as: 


Tn+1  Tn  )  aAt 

Yj  '  Ti.j  exp  1 77772 


(4x) 


Tn  Tn  _n  Tn  .,n 

TU1J  T1-U  *  Ti,j+i  Ti  J-i  "  4TU 


r-n 


1.J 


-7(11) 


Note  the  At  that  appears  In  the  above  Is  the  time  elapsed  between 


time  steps  n  and  n+1 .  The  temperatures  on  the  right-hand  side  are 

>  th 


the  four  nearest  neighbors  to  the  1,j  node  (see  Fig.  1). 

In  keeping  with  the  notation  derived  by  Bhattacharya  [2],  the  term 

aAt 


«  = 


(  AX) 1 


(12) 


Is  called  the  dimensionless  time  step  and 


1 .  j 


Tn  Tn  Tn  Tn 

1+1*1  *  'l-l.J  *  *  'l.J-1 


4T 


LI 


(13) 


1.J 


Is  called  the  dimensionless  drive  number.  Thus,  Eq.  (11)  can  be 
rewritten  rather  compactly  as: 

,n+l  , n 


{"?.!>  (,4) 

It  was  found  [2]  that  an  Improvement  In  the  solution  at  the  n*l 
time  step  at  node  (1,J)  can  be  made  If  the  time  step  Is  divided  Into 
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a  number  of  equal  length  time  sub-intervals.  The  method  of  time  step 

sub-intervals  can  be  described  as  follows.  Let  us  assume  that  the  time 

Interval  will  be  divided  Into  three  Intervals  Including  the  one  at  the 

end  of  the  time  step  (Fig.  2).  Now  returning  to  Eq.  (14),  and 

evaluating  at  the  n+1/3  time  step,  results  In: 

n+1/3  Tn  a  IG..n  \  nM 

TUJ  sTu  "'ts'ul  (15) 

Proceeding  from  the  n+1/3  time  step  to  the  n+2/3  time  step: 

Tn+2/3  Tntl/3  „,J5ynt1/3l  ,  i  v 

T1J  ’  TU  eXP\3Ht.J  /  (,6) 

Finally,  the  temperature  can  be  written  at  the  n+1  time  step: 

,n+l  rn+2/3  .wJ2yn+2/3\  mi 

T1.J  =  T1.J  exp\3MAJ  ) 

Now  substituting  Eq.  (15)  Into  Eq.  (16)  and  substitution  of  Eq.  (16) 
Into  Eq.  (17),  produces: 

Tn+1 


1.3 


r  a  <  Q  Mn  \  /£  „n+l/31  i  a  Mn+2/3\ 

f1.J  CXP\3  M1,j}  exp \3  M1J  /eXP\3M1,J  / 


(18) 


or 


+  m 


n+1/3 

1.3 


(19) 


where  the  M's  are  then  evaluated  at  the  sub-time  Intervals  and  then 
summed  for  calculation  of  "T"  at  the  n+1  time  step.  Equation  (19) 
can  be  written  more  generally  as: 


Tn 

Ti,j  exp 


m 

ft  \  '  un+p/(m+l) 
m  +  1  Lt  1.3 

p=0 


(20) 


In  reference  [2],  It  was  shown  that  for  heat  transfer  applications,  the 
time  step  can  be  subdivided  as  follows: 
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j  ?f  -  1  .  .  .  heat  transfer  coefficient  -*  ® 


(21) 


[  |  +  1  .  .  .  finite  heat  transfer  coefficient 
’  1  with  flanking  nodes  outside 

calculation  domain. 

N  =  number  of  nodes  In  one  of  the  coordinate  directions. 

The  procedure  necessary  to  determine  the  dimensionless  drive 
numbers  will  now  be  described.  Since  the  method  Is  an  explicit 
technique,  all  information  Is  known  from  the  previous  time  step  or  from 
the  previous  time  sub-interval.  The  effort  Is  the  centered  around 
calculation  of  the  drive  numbers  at  the  requested  number  of  time 
sub-intervals  for  each  of  the  spatial  positions  (nodes). 

The  calculation  procedure  Is  shown  In  Fig.  (3).  Because  the  drive 
numbers  are  evaluated  at  sub-time  Intervals,  the  temperature  (or  any 
other  field  variable)  must  also  be  known  at  these  sub-time  Intervals. 
Therefore,  the  temperature  field  Is  calculated  at  each  sub-time 
Interval,  and  In  turn  Is  used  to  calculate  the  next  drive  number.  The 
drive  numbers  for  each  node  are  summed  for  all  the  sub-time  Interval 
steps  and  then  used  to  predict  the  field  variable  at  the  next  complete 
time  step.  This  results  In  a  computer  storage  requirement  of  4(N), 
where  N  Is  the  number  of  nodes. 


Extension  To  One-Dimensional  Cylindrical  Coordinates 
In  the  previous  section,  the  exponential  finite  difference 
technique  was  extended  to  two-dimensions.  Now  the  procedure  will  be 
considered  In  another  coordinate  system.  In  particular  the  method  will 
next  be  applied  to  radial  one-dimensional,  unsteady  heat  transfer. 

The  governing  unsteady  diffusion  equation  for  a  material  with 
constant  properties  Is: 
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(22) 


'Cp  S  -  k 


[iL 

(r 

[r  a r 

\  ar/J 

alt  + 

1  111 

ar2 

r  ar 

Assume  that  the  temperature  can  be  written  as  the  product: 

T(r  ,t)  =  «|»(  r  )o(  t) 

The  Initial  conditions  are 

T(r,0)  =  f(r)  ;  9(0)  =  1 

Now  taking  the  appropriate  derivatives  of  Eq.  (23);  results  In 


(23) 


(24) 


ai  ae  .  aj  A  a*  .  lit  A 


at 


at  •  ar 


ar  * 


(25) 


ar 


ar 


Substituting  Eq.  (25)  Into  Eq.  (23)  yields: 

.2. 


4» 


'  2  T 

as  A  a  9  1  a» 

at  =  a  [e  *  F  0  a? 


Olvldlng  both  sides  by  90  brings: 


1  ao  _  1  a29  l  99 

e  at  =  a  9  2  r9  dr 

di 


9  2 
w  ar 


(26) 


It  can  be  seen  that  the  variables  have  been  separated.  Thus,  both 
sides  of  Eq.  (26)  must  equal  a  constant,  k.  Now  using  the  time 
dependent  side  of  the  above  ,  l.e.. 


1 

0  at 


and  multiplying  this  by  09/09  produces: 


1  ai 

T  at 


-  tc 


Integrating  for  a  particular  value  of  the  radial  position 

T  *  C]  exp  {-  ict } 


"r"  gives: 


*  -  «  #.<  ii*  !»*  #.»•>,<  «.*  >.* 
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Next  the  Initial  condition  Is  used  and  the  equation  can  be  written  as: 

T(r.t)  =  T(r,0)  exp  {-  *t}  (27) 

Returning  to  Eq.  (26)  and  using  the  radial  side  of  the  equation,  we  have 


1  3  $  +  1_  3£  1 

|_*ar2  r*arJ 


-  K 


ar2  r* 

Multiplying  by  64>/e<t>,  the  equation  can  be  rewritten  as: 


|~1  a2T  1_  ill 

LT  ar2  TrarJ 


(28) 


Incorporating  the  Initial  condition, 


T( 


o _ fa2!  ,  1  aTl 

r’° )  Lar2  r  arJ  = 


(29) 


Next  the  partial  derivatives  are  replaced  using  central  finite 
differences  [4]: 


HI 


+  T 


1-1 


2T 


(Ar)1 


aj 

ar 


HI 


2  Ar 


n 

1-1 


Substituting  Eq.  (30)  Into  Eq.  (29): 


(30) 


(31) 


Equation  (31)  Is  used  to  replace  the  constant,  -  k  In  Eq.  (27),  l.e.. 


rn+l 


t"  exp 


At 


Tn  *■  T° 
'HI  1 1-1 


-  2T 


(Ar)' 


yn  _  Tn 

'HI  1 1-1 

2  Ar 


Rearranging  this  brings 


Equation  (32)  states  that  the  temperature  at  the  1  radial 
position  at  the  (n+1)  time  step  Is  found  from  the  product  of 
temperature  at  1th  position,  n*h  time  step  and  the  exponential 
term  that  Is  composed  of  a  dimensionless  time  step: 


and  a  dimensionless  drive  number: 


It  should  be  noted  that  this  drive  number  applies  to  the  Interior  nodes 
(r^  *  0).  For  the  node  at  r  »  0  the  dimensionless  drive  number 


becomes: 


»W.,  -  0 


Finally  Eq.  (32)  can  be  written  In  a  compact  form  as 

t"+1  =  exp  }  (36) 

The  sub-time  Interval  evaluation  for  Eq.  (36)  Is  the  same  as  that  found 
In  the  two-dimensional  Cartesian  form  as  shown  earlier. 


Stability  of  the  Exponential  Finite  Difference  Method 


With  few  exceptions,  explicit  finite  difference  procedures  for 
solving  partial  differential  equations  are  Inherently  unstable,  unless 
certain  numerical  condltons  are  satisfied.  These  conditions  take  the 


w 


y. 


form  of  a  grid  size  and/or  time  step  requirement  written  In  terms  of 
parameters  of  the  given  problem.  If  these  stability  conditions  are  not 
met,  the  solution  will  diverge,  often  rather  drastically.  On  the  other 
hand,  the  stability  requirements  can  make  explicit  methods  Impractical 
for  a  particular  application  by  requiring  an  unrealistically  small  grid 
or  time  step.  Nonetheless,  these  conditions  must  be  known  prior  to  the 
use  of  any  explicit  finite  difference  procedure. 

There  are  a  variety  of  methods  that  have  been  used  to  establish 
the  stability  constraints  of  a  finite  difference  procedure:  some  are 
very  elementary,  some  quite  Involved.  In  essence,  the  methods  seek  to 
find  an  expression  for  the  amplification  factor  which  Is  the  ratio  of 
the  current  solution  result  to  that  In  the  previous  step.  If  the 
absolute  value  of  the  ratio  Is  less  than  one  the  method  Is  stable. 
Determination  of  the  amplification  factor  for  the  exponential  finite 
difference  method  Is  particularly  convenient,  as  has  been  shown  In 
[2].  For  the  one-dimensional  cylindrical  coordinate  case,  the 
amplification  factor  i  can  be  readily  defined  as 


or  from  Eq.  (32)  for  an  Interior  node: 


So  for  stability  to  exist  as  at  and  Ar  approach  zero: 
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11m  m  <  1  (39) 

At  -»  0 
Ar  -*  0 

To  satisfy  this  requirement  the  exponent  In  the  eponentlal  of 
Eq.  (38)  must  obviously  be  less  than  or  equal  to  zero.  Since  the 
components  that  make  up  ft  In  that  exponent  are  all  positive,  this 
Implies  that  the  dimensionless  drive  number  will  dictate  whether  or  not 
the  stability  criteria  Is  met.  For  the  cylindrical  coordinate  case  the 
dimensionless  drive  number  must  satisfy: 

Multiplying  by  t”  this  becomes: 

(TU.  *  T?-l  -  2T?)  •  TT,  (TU!  -  T?-l)  i  ° 

Define  B  =  Ar/2r^  and  rearrange  to  get 

(1  ♦  0)T"+1  ♦  (1  -  0)Tj  <  2tJ 

(41) 

T"  >  \  (1  +  B)T"+1  *  (1  -  B)T"  n 

Equation  (41)  needs  to  be  satisfied  otherwise  an  unstable 
condition  can  exist.  In  Eq.  (41)  as  Ar  -*  0,  or  equivalently  0  -»  0, 
the  stability  condition  becomes: 

Tn  >  1  -rn  ♦  Tn  (42) 

1-2  'm  '  1-1 

Also  for  one-dimensional  cylindrical  coordinate  case,  the  node  at 
r  ■  0  has  a  different  stability  criteria  because  the  node  at  1*1 
does  not  exist.  Since  the  radial  derivative  of  the  field  variable  must 
equal  zero  at  r  =  0.  In  finite  difference  terms  this  can  be  realized 


a'V'.'s’ 


»*r  t'lj'lj  «.l  .  < 


by  requiring  that  »  T^  -j  at  this  node.  Thus,  at  the 

origin,  Eq.  (42)  becomes 


Tn  >  7^ 

'l  -  1 1-1 


As  stated  by  Bhattacharya  [2],  the  dimensionless  drive  number  Is 
the  determining  factor  whether  or  not  the  stability  criteria  Is  met. 
However,  the  dimensionless  time,  If  made  large  enough,  could  cause  the 
solution  to  become  unstable.  Since  time  sub-interval  division  Is  used, 
the  total  dimensionless  time  step  Q  could  become  quite  large.  From 
[2]  It  was  recommended  that  the  dimensionless  time  step  be  chosen  to 


satisfy  the  following: 


m  ♦  1 


<  0.5 


where  m  Is  the  number  of  time  step  sub-intervals  Involved  In  the 
calculations.  If  Q  *  5,  for  example,  then  m  would  have  to  be 
greater  than  or  equal  to  9  l.e.,  nine  sub-intervals  would  be  required. 
From  Eq.  (21)  with  Infinite  heat  transfer  coefflcent,  a  minimum  of  20 
nodes  would  be  needed. 

Another  useful  comparison  to  be  made  Is  the  one-node  model  as  used 
In  Refs.  [2]  and  [7]  where  the  value  of  the  dependent  variable  at  the 
surrounding  nodes  Is  set  equal  to  zero.  For  the  one-node  model,  as 
stated  In  Ref.  [7],  the  exponential  finite  difference  and  exact  are  the 
same  (Fig.  4).  This  figure  Indicates  that  the  exponential  solution 
remains  stable  as  Q  Is  Increased. 
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Effect  of  Initial  and  Boundary  Conditions  on  the  Exponential 
Finite  Difference  Method 

In  this  section  the  effect  of  boundary  conditions  on  the 
exponential  finite  difference  technique  will  be  Investigated.  Boundary 
conditions  that  are  typical  of  heat  transfer  applications  will  be 
considered.  The  conditions  to  be  presented  are:  (1)  finite  heat 
transfer  coefficient  (mixed  condition),  (2)  Infinite  heat  transfer 
coefficient  (Dlrlchlet  condition),  (3)  constant  heat  flux  (Neumann 
condition),  and  (4)  time  varying.  Initial  conditions,  where  the  field 
variable  Is  equal  to  zero,  will  also  be  discussed. 


Finite  Heat  Transfer  Coefficient 

For  this  boundary  condition,  the  method  used  In  reference  [2]  will 
be  utilized.  Numerical  Implementation  of  the  boundary  condition 
requires  that  a  node  be  placed  outstde  the  solid  In  the  surrounding 
medium.  This  external  node  will  be  used  In  the  finite  difference 
equation  at  the  solid  surface.  One-dimensional  heat  transfer  In  a 
cylinder  (T  -  T(r,t))  will  be  used  to  demonstrate  the  procedure. 

Using  the  exponential  finite  difference  method  for  T  -  (r,t).  It 
was  found  earlier  (Eq.  (36))  that 

T^1  =  t"  exp  |«mJJ  (45) 


•w 
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where  the  dimensionless  time  Is  given  by 


ft 


a  At 
(&r)2 


and  the  dimensionless  drive  number  by: 

rn' 


_n  Tn 
TU1  +  T1-l 


2T 


1 


(46) 


The  thermal  condition  at  the  surface  Is  found  by  equating  the 
conductive  and  convective  heat  fluxes  at  the  surface: 


(47) 


Using  the  node  numbering  as  shown  In  Fig.  5  and  a  central  difference  to 
express  the  deveratlve,  Eq.  (47)  becomes: 


Solving  for  t|J,  the  temperature  of  the  external  node,  yields 

,n  ,n  2h  dr  /_  ,n 
T0  -  T2  *  ~  (T-  *  T1 


(48) 


(49) 


Let  B  -  h  dr/k  (Blot  number  [8]),  thus  Eq.  (49)  Is  written: 

Tjj  x  ♦  2BT^  -  2BlJ  (50) 

Now  that  an  expression  for  the  temperature  at  the  external  node  Is 
known.  It  can  be  used  In  the  expression  for  the  temperature  at  the 
surface  node.  This  results  In: 


Tn+1  Tn  „  /_  Mn) 

T1  =  T1  exp  (ft  M1  J 

where  (51) 


_n  __n  _n  „„_n 

T«  -  2T.  ♦  2BT  *  T 0  -  2BT. 

2  1  ®  2  1 

[«(t:  -  t?) 

J 

♦ 

+  R 

tJ 

L  l 

i ^  if-  fiV.'.  A  A  f. 


'■rlylvLv.v. 


l*.V k  kZjC x* Ja  x  *  ■ 


and  the  drive  number  can  be  further  simplified  to: 


M 


n 

1 


2  t[!  -  (1  ♦  B)Tj  ♦  BT.l 

- 1 

c  »— 
H— 

_ 1 

ArB 

R 


(52) 


Infinite  heat  transfer  coefficient 
If  the  heat  transfer  coefficient  In  Eq.  (47)  Is  placed  on  the 
left-hand  side  of  the  equation  and  then  allowed  to  become  Infinite,  It 
Is  seen  that  the  surface  temperature  will  equal  the  temperature  of  the 
surroundings.  Thus,  the  boundary  condition  In  which  h  -»  ®  Is 
Identical  to  that  In  which  a  boundary  temperature  Is  held  constant.  In 
the  calculation  procedure,  these  Isothermal  boundary  nodes  are  only 
needed  for  calculation  of  the  temperature  field  at  the  surrounding 
nodes.  For  example  In  a  two-dimensional  square  grid  with  a  total  of 
121  nodes  calculation  would  be  reduced  to  a  total  of  81  nodes  If  the 
temperature  Is  specified  for  all  four  boundaries. 


Constant  Heat  Flux 

For  a  constant  heat  flux  applied  to  the  boundary  surface,  the  same 
procedure  as  was  used  In  the  finite  heat  transfer  coefficient  case  will 
be  utilized.  The  condition  at  the  surface  for  one-dimensional 
cylindrical  coordinates  Is  given  by: 

An  external  node  Is  placed  outside  the  solid  as  was  done  earlier. 
Equation  (53)  can  then  be  written 
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as  the  subtime  Interval  calculations  are  made. 

Dependent  Variable  Initially  Equal  to  Zero 
One  last  condition  that  can  exist  wherever  there  Is  Initial  zero 
temperature  (T(xf0)  =  0)  needs  to  be  discussed.  If  this  condition  Is 
encountered,  then  the  following  substitution  should  be  made  or  else  the 
exponential  finite  difference  method  will  not  work.  This  can  be 
readily  seen  by  examlng  any  of  the  numerical  equations  e.g.,  Eq.  (56). 
Since  the  Initial  temperature  would  appear  In  the  denominator  of  the 


—— a— — — — MTBWT— 1 BBHBMWagnui  Hi  WJ  PUI  HUM 'Ml 
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exponent  In  the  exponential,  problems  would  ensue.  To  circumvent  this 
difficulty  define  a  new  variable  T,  such  that  T(x,t)  >1.0-  T(x,t). 
Now  the  exponential  finite  difference  equations  described  above  can  be 
utilized  by  simple  replacement  of  the  1  variable  with  the  T 
variable. 


a 

g 


*-*'*.*  JAM  r*  f*V  f  ■>  1 1 '  It  .1*  .«|  a  I  a’l  «  I  m  l.l'l.l 
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III.  NUMERICAL  COMPARISON  OF  THE  EXPONENTIAL  FINITE  DIFFERENCE 
METHOD  TO  EXACT  SOLUTIONS  AND  OTHER  NUMERICAL  TECHNIQUES 
The  final  product  of  any  numerical  study  Is  how  well  the  given 
method  performs  when  compared  to  known  exact  solutions  or  to  other 
numerical  techniques.  The  exponential  finite  difference  method  will 
now  be  applied  to  the  following  cases  to  demonstrate  the  capability  of 
the  method  to  solve  the  diffusion  equation: 

(I)  One-dimensional  heat  conduction  In  cylindrical  coordinates  with  an 
Infinite  and  a  finite  heat  transfer  coefficient  at  the  surface, 
unsteady  state, 

(II)  Two-dimensional  unsteady  state  heat  conduction  In  Cartesian 
coordinates, 

(III)  Solution  of  Laplaces  equation,  Cartesian  coordinates, 

(1v)  One-dimensional  heat  conduction  In  Cartesian  coordinates,  unsteady 
state  with  temperature  varying  thermal  conductivity, 

(v)  Steady  state  Couette  flow, 

(vl)  Three-dimensional  heat  conduction  In  Cartesian  coordinates, 
unsteady  state. 


One-Dimensional  Heat  Conduction  In  Cylindrical  Coordinates 
The  one-dimensional  cylindrical  coordinate  heat  conduction  case 
with  temperature  as  a  function  of  time  and  radial  position  will  be 
Investigated  for  Infinite  and  finite  heat  transfer  coefficient.  The 
exact  results  for  both  cases  can  be  found  In  Ref.  [9]. 
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For  Infinite  heat  transfer  coefficient  on  the  boundary  surface  the 
exact  result  Is  given  In  [9]  as: 


T(r.t)  -  T 


=  2 


-aX?t 

e 

m=l  (V)  Jl<XmR> 


Z 


(57) 


where  X-R  Is  the  zero  of 

J0(X-R)  =0  (58) 

The  results  of  both  the  exact  analysis  and  the  exponential  finite 
difference  method  are  shown  In  Table  I.  As  can  be  seen  from  the 
tabulated  results,  exponential  finite  difference  results  approach  the 
exact  solution  as  the  number  of  nodes  Is  Increased  or  as  the 
dimensionless  time  step  Is  decreased. 

When  the  heat  transfer  coefficient  has  a  finite  value  at  the 
surface,  the  exact  solution  from  [9]  Is: 


T(r,t)  -  T. 


2B 


(59) 


where  B  =  hR/k  (Blot  number)  and  X-  (characteristics  values)  are 
given  by  (for  cooling): 

<V> W>  -  BVV>  ■ 0  (60) 

The  results  are  shown  In  Table  II  for  various  values  of  the  Blot 
number.  As  would  be  expected  the  solution  approaches  the  exact 
solution  as  the  number  of  nodes  In  Increased.  The  size  of  the  Blot 
number  did  not  seem  to  effect  the  accuracy  of  the  solution.  As  the 
elapsed  time  of  the  solution  proceeded,  temperatures  predicted  by  the 


exponential  finite  difference  method  approached  the  exact  result.  Also 
the  results  Indicated  that  reducing  the  size  of  the  time  sub-interval 
Increased  the  method's  accuracy. 

One  last  comparison  will  be  made  while  Investigating  the 
exponential  finite  difference  technique  In  one-dimensional  cylindrical 
coordinates.  The  problem  situation  Is  shown  In  Fig.  (6)  and  applies  to 
a  cylindrical  annulus  with  the  following  Initial  and  boundary 
conditions : 

T( r ,0)  =  0 

T(R2.t)  =  1.0  (61) 

If  <"!•*)  •  0 

In  Ref.  [4]  this  problem  was  solved  numerically  using  a 
characteristic-value  solution.  A  comparison  of  results  Is  shown  In 
Table  III  for  the  exponential  method  using  the  same  grid  spacing  as  In 
[4]  and  for  the  case  where  grid  spacing  Is  halved.  The  results  are 
seen  to  compare  quite  well  with  the  finer  mesh  being  slightly  closer  to 
the  value  from  Ref.  [4]  especially  during  the  first  few  time  steps  of 
the  solution. 

Two-Dimensional  Heat  Conduction  In  Cartesian  Coordinates 

The  exponential  finite  difference  technique  will  now  be  applied  to 
the  two-dimensional  heat  conduction  problem  In  Cartesian  coordinates 
shown  In  Fig.  (7).  The  exponential  finite  difference  method  will  be 
compared  to  the  solution  of  this  problem,  as  performed  In  Ref.  [4], 
using  the  alternating  direction  Implicit  technique  (ADI).  The  results 
of  the  two  numerical  techniques  and  the  exact  analysis  are  shown  In 
Table  IV.  The  temperature  Indicated  for  comparison  Is  that  at  the 


origin  x  *  y  =  0,  shown  In  Fig.  (7).  As  maybe  seen,  the  ADI  technique 
does  not  predict  the  temperature  as  accurately  as  the  exponential 
finite  difference  method  at  the  first  time  value  shown  In  Table  IV. 
However,  as  the  amount  of  elapsed  time  Increases  either  method  does  a 
very  good  Job  at  predicting  this  temperature.  When  the  number  of  grid 
points  was  Increased,  by  halving  the  spatial  Intervals,  the  exponential 
finite  difference  method  was  found  to  be  more  accurate  for  all  the  time 
steps. 

Since  the  ADI  method  Is  one  that  requires  simultaneous  solution  of 
equations  In  the  two  coordinate  directions,  the  time  step  size  can  be 
made  large.  The  exponential  finite  difference  technique  must  have  the 
dimensionless  time  step  kept  below  0.25  to  keep  the  solution  stable. 

So  the  required  CPU  time  for  the  exponential  method  Is  higher  for  this 
application. 

Solution  of  Laplaces  Equation 

Since  the  exponential  finite  difference  method  has  been  used  for 
two-dimensional  unsteady  state  conduction,  a  natural  extension  with 
little  additional  effort  would  be  to  use  this  method  to  solve  Laplace's 
equation.  This  can  be  Implemented  In  the  exponential  finite  difference 
method  by  Just  allowing  the  solution  to  march  In  time  until  no  further 
change  In  the  field  variable  Is  Indicated. 

As  an  example,  the  problem  as  shown  In  Fig.  (8)  will  be  solved  and 
the  results  compared  to  those  given  In  Ref.  [4J.  In  the  referenced 
work,  the  solution  was  found  by  using  a  Gauss-Seldel  Iterative 
technique. 


A  comparison  of  results  along  a  diagonal  from  the  position  (x  =  0, 
y  >  1)  to  (x  >  1,  y  »  0)  Is  presented  In  Table  V  for  two  different  grid 
spaclngs.  As  can  be  seen,  the  solutions  are  nearly  Identical  with  the 
exponential  method  requiring  a  smaller  number  of  Iterations  (or  time 
step  Increments)  to  reach  a  slmlllar  result. 

One-Dimensional  Unsteady  State  Conduction  With  Temperature 
Varlng  Thermal  Conductivity 


The  effect  of  temperature  varlng  thermal  conductivity  will  now  be 
Investigated  using  three  different  numerical  schemes:  a  pure  explicit, 
the  exponential  method  and  an  Implicit  technique.  The  problem  to  be 
solved  Is  Illustrated  In  Fig.  (9a).  The  thermal  conductivity  as  shown 
In  Fig.  (9b)  Is  assumed  to  be  a  linear  function  of  temperature. 

The  exponential  finite  difference  method  will  be  applied  first  to 
the  given  problem.  The  governing  partial  differential  equation  Is  [1]: 


“p  at  '  ax 


(*S) 


From  [1],  Eq.  (62)  can  be  changed  to  a  simpler  form  by  using  a  new 


variable  e  (the  Klrchoff  transformation)  given  by: 


k(T)dT 


(63) 


where  kR  Is  the  conductivity  at  temeprature  TR,  and 

ae  k_  aj  aj  _R  ae 

at  "  kD  at  or  at  =  k  at 

R  (64) 

ae  _  k_  al 
ax  =  kR  ax 


Substituting  Eq.  (64)  Into  (62)  gives: 
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/ae\  a__  (u  ae\ 
\at/  *  ax  yKR  ax ) 


oC  /  \  2 

_fi  /a©\  a_e 

k  UtJ=ax? 


Since  It  has  been  assumed  that  the  thermal  conductivity  Is  a  linear 


function  of  temperature. 


k(T>  =  kR(l  *  (IT) 


Now  returning  to  Eq.  (63)  and  substituting  In  Eq.  (66),  we  have: 


kl  <k*‘ 


BkRT)dT 


Direct  Integration  yields: 


9  -  <T  -  v  {'  ♦  I  <T  *  V} 


Equation  (67)  provides  the  relationship  between  the  variable  T  and 


the  new  variable  e. 


Returning  to  Eq.  (65) , 


ae  k  a  e 


3t  '  O'p  Sx2 


Equation  (68)  Is  In  a  form  now  that  the  exponential  finite 


difference  can  be  applied.  The  resulting  equation  In  the  variable  can 


be  shown  to  be  given  by: 


n.l  n  at  kt  (eUl  *  91-1  ‘  29l) 

■ e' exp  ^  [ — >  ; 


tm 
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Evaluating  Eq.  (67)  at  node  1  brings 


8 ?  ■  K  -  tr)  I'  *  I  (Tt  *  tr ) j 

•  m  -M  *i  (M)2  *  i) 


Substitution  of  Eq.  (70)  Into  Eq.  (69)  at  the  appropriate  time  steps 


and  nodal  locations  will  give: 


exp  jyk1 


where 


M*!(K*’)2-r)=  [(’M«H[K)2- 
n  [«.,  *  MMR/  *  MJMHF 


2  - 12  . 

R 


M  -  tr)  *  I  [ft) 


2  -  T2 

R 


pCp(&x)‘ 


If  Tr  =  T,,,  =  0.0,  Eq.  (71)  becomes: 


-r1  * !  ft*’)2  ■  (M  ft)  • 

I . n[K.i :  ftJ- ft  *  f  [KJ2  *  ftJ2  -  2ft): 

r  l  ft*  i«)2 


The  equation  for  T^  Is  a  quadratic  with  the  right-hand  side  of  the 
equation  all  being  known  at  time  step  n,  so  define  a  variable  ^ 
such  that 

•t  •  [ft  ♦ !  (ft)2] ex|1  ft* 

(fto : ft.i  - ft) .  I[(t"m)2  *  (ftft2  -  2(ft)2]]  m. 

[  ’MM)2  J 


27 


Equation  (72)  then  becomes: 

/Tn+1  \2  2  Tn+1  2  n 
lT1  >  *  8  Tt  -  8  ”t  *  0 

Solving  this  and  using  the  positive  root  results  In: 


(74) 


r  •  s  (-  1  *  V  *  2«,«  ) 


(75) 


where 

B  >  0 

Equation  (75)  In  conjunction  with  Eq.  (74)  are  Implemented  In  the 
exponential  finite  difference  solution  sequence.  In  this  case  the 
conductivity  as  well  as  the  temperature  field  must  be  kept  track  of  on 
the  sub-time  Interval  level.  The  dimensionless  time  step,  Q,  and  the 
rate  of  conductivity  change,  B,  must  be  both  considered  when  choosing 
the  step  size  so  the  solution  does  not  become  unstable.  For  this 
method,  the  term  (yklj/m+l)  in  the  exponential  was  considered  at 
Its  maximum  possible  value  and  the  time  step  was  adjusted  to  retain 
stability.  This  criteria  was  chosen  so  that 

n 


Y  k 


1 


<  0.5 


m  +  1 

The  next  method  to  be  Investigated  for  the  temperature  varying 
conductivity  problem  will  be  the  pure  explicit  method.  As  stated 
earlier  the  governing  partial  differential  equation  for  one-dimensional 
conduction  Is  given  by: 

3T  a 


pc. 


(*H) 


' p  at  ~  ax 

Using  the  chain  rule  this  equation  can  be  put  Into  a  nonconservative 


(76) 


form. 


ESQ 
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.  al  .  t  8k  /dl\2 

p  P  at  *  *  2  aT  [ax } 

'  ax  \  / 


Assuming  the  same  linear  profile  as  before: 

k(T)  =  kR  ♦  QkRT 


—  -  Bk 
3T  "  WR 


Substituting  Eq.  (78)  Into  (77)  will  give: 


'CD  S  ■  V1  ‘  0T> 


ill  ,  Bk  /iiv 

ax2  P  R  \3x/ 


S  •  v1  *  M>  £ 


1  B  /ai\2 

2  BaR\axj 


where 


R 

Using  central  space  differences  and  a  forward  time  difference  we  may 


write: 


jn  _  yn 

3T  'HI  '  1-1 


Tn+1  Tn 
aT  *1  ~  *1 


2  Tn  +  Tn  -  2Tn 
3  T  'HI  1 1-1  *'l 


Substituting  Eq.  (79)  Into  (78)  produces: 
n+1  n  [Tn  n  Tn  1  fTn 

- it - “si'  *  BT1  I  ^2  0  "  l 


-  T 

+i _ Vj 


S'fVJyV 


y/yy-y  «> y  y  y  y  .  y : 


I 


This  may  be  simplified  and  written  as: 


rn+l  xn  ~  IT.  OTn  \  /Tn  _n  _Tn\  |  /Tn  Tn  \2* 

F1  "  T1  +  0  |^\1+BT1  )(T1+1  +  T1-l  ~  2Tl)  +  4  (T1+l  "  Vlj 


where 


aR  At 


Equation  (80)  can  be  used  to  directly  solve  for  the  temperature 


field  at  the  next  time  step.  Some  additional  care  must  be  used  to  keep 


the  solution  stable  as  the  size  of  the  dimensionless  time  step  Q, 
and  the  rate  of  conductivity  change,  0  both  effect  the  solution. 

The  final  method  to  be  Implemented  for  comparative  purposes  Is  the 
Implicit  method.  Starting  with  Eq.  (78),  we  have: 


ft  ■  *R  ('  *  BT>  *  %  (0 


To  avoid  a  solution  sequence  that  would  require  the  solution  of 


nonlinear  algerbralc  equations,  the  following  will  be  assumed; 

(a)  The  term  (l  *  0T)  can  be  replaced  with  (1  *  Bt”) 

(b)  The  squared  first  derivative  term  can  be  replaced  by 


/  \2  An  -  Tn  \  /t"+1  -  Tn+1 

(21)  1  +  1  Vj  (ll±l  1-1 

\dx/  \  2Ax  I  \  2Ax 


This  Is  a  linearizing  technique  known  as  lagging  the  coefficients. 
Substituting  these  Into  Eq.  (81)  will  produce: 


rn+l  Tn 


—  *  i  -n 

Ti  ~  T<  L 
at  "  “r  \ 


1  +  0T 


Tn+1  n+1  9Tn+l 

Jl+1  +  1-1  ~  1 


+  aR0 


n  Tn  \  /rn+1  Tn+1 

1+1  '  1 1-1  )  [ 1 1+1  "  1 1-1 


Further  simplification  gives: 


;rtrw.r 


w.v.\;v  mv  .vvv.v.v. 


£ 


j] 


Tn+1  Tn 
'l  “  'l 


where 


Now  define: 


»(’  *ST?)(T? 


n+1  Tn+1 
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Substituting  Eq.  (83)  Into  (82)  yields: 


Tn+1  Tn 
'l  "  'l 
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Simplifying  this  results  In: 
T 


^1  +  2Qn”  Y  (onj  ♦  flkjj  -  t£{  -  CkJ  j  -  tJ  (84) 

The  equation  shown  above  Is  now  In  a  form  that  can  be  used  In  the 
Thomas  Algorithm  [4],  l.e.,  Eq.  (84),  can  be  written  as: 

.n 


where 


_n+l  .  Tn+1  _n+l 
a1T1-l  *  V)  +  c1T1+l 


a>  ■  -  °K  - x?) 

b1  =  1  +  2«n" 

c,  -  -  a(»J  *  x?) 
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(85) 


Equation  (85)  can  now  be  solved  using  a  tri-diagonal  matrix  routine. 
The  variables  &y  b^,  and  c^  must  be  evaluated  at  each 
position  and  time  step  as  their  values  change  as  the  field  variable  T 
changes. 


A  comparison  of  results  of  the  three  methods  can  be  found  In 
Fig.  (10)  and  Table  VI.  Figure  (10)  shows  the  temperature  field 
through  the  slab  cross-section.  From  this.  It  Is  evident  that  the 


exponential  and  pure  explicit  methods  give  very  slmlllar  results.  The 
Implicit  method  predicted  higher  temperatures  closer  to  the  slab 
surface  and  lower  temperature  at  the  slab  centerline  then  either  of  the 
two  explicit  methods.  In  Table  VI  the  results  at  the  slab  center  are 
shown  for  various  elapsed  times.  As  can  be  seen,  all  three  methods 
agreed  with  each  other  to  within  a  few  percent. 


The  Steady  State  Temperature  Field  of  a  Couette  Flow 
Another  application  of  the  exponential  finite  difference  method 
will  now  be  presented.  The  problem  to  be  Investigated  Is  the 
developing  temperature  field  In  laminar  Couette  flow  [7].  The  problem 
statement  Is  illustrated  In  Fig.  (11).  Neglecting  viscous  effects,  the 
governing  equation  Is  given  by  [10]: 


pC 


U  —  = 

p  x  ax 


(86) 


Neglecting  conduction  In  the  x-dlrectlon  or  assuming  that  the 
convection  term  Is  much  greater  than  the  conduction  term,  Eq.  (86) 
becomes : 


U 


31 
X  dX 


a 


From  Fig.  (11)  using  the  expression  for  the  velocity  In  the 
x-dlrectlon,  Ux  =  Uy/L,  Eq.  (87)  becomes 


(87) 
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ai  .  a_I  (88) 

ax  Uy  ay2  v  ; 

Equation  (88)  Is  now  In  a  form  where  separation  of  variables  can 
be  Implemented.  Following  the  same  procedure  as  Indicated  earlier  to 
find  the  exponential  finite  difference  equation.  It  can  be  shown  that: 

Tj-  .  Tj  eKP  !li\  <89, 

0  J  U(&y  Tj  / 

The  procedure  utilized  here  Is  that  the  solution  marches  In  the 
x-dlrectlon  Instead  of  time  as  was  the  case  for  the  previous  examples. 
Information  from  the  last  x-posltlon  and  y-dlrectlon  are  used  to 
determine  the  dependent  variable  at  the  next  x-posltlon. 

Results  of  Implementing  this  method  are  shown  In  Fig.  (12).  The 
temperature  field  Is  shown  for  three  x-locatlons  for  two  different 
values  of  L/U.  The  results  Indicate  that  as  the  upper  plate  velocity 
U  Is  Increased,  the  propogatlon  of  the  temperature  change  In  the 
y-dlrectlon  Is  slowed  down. 


Unsteady  State  Heat  Conduction  In  Three-Dimensional  Coordinates 


The  final  application  of  the  exponential  finite  difference  method 
to  the  diffusion  equation  will  be  that  of  three-dimensional,  unsteady 
state  heat  conduction.  The  exponential  method,  a  pure  explicit  method, 
and  an  Implicit  method  (method  of  Oouglas,  [11])  will  be  compared  to 
exact  solution  for  the  situation  shown  In  Fig.  (13). 

The  exact  solution  to  the  problem  Illustrated  In  Fig.  13  Is  given 
In  Ref.  [10]  as: 


T(x,y,z,t) 


•  cos [(”  *  2) fjcos[(n  *  2) “]cos[(p  *  2)?]  (90) 

where  a,  b,  and  c  are  the  widths  of  the  cube  In  the  x,  y,  and  z 
directions  respectively.  Equation  (90)  will  be  used  to  determine  how 
well  the  numerical  techniques  predict  the  temperature  distribution. 

The  exponential  finite  difference  technique  will  be  Investigated 
first.  The  sequence  to  be  followed  for  determining  the  finite 
difference  equation  Is  the  same  as  presented  for  the  earlier  cases. 

The  procedure  for  this  three-dimensional  case  consists  of  the  following 
stepped  procedure: 

(1)  Linearize  the  partial  differential  equation 

(2)  Assume  a  product  solution 

(3)  Separate  time  from  spatial  dependence 

(4)  Solve  for  time  dependence 

(5)  Insert  the  appropriate  spatial  finite  differences  Into 
exponential  term  that  results  from  step  3 

Based  on  this  procedure  the  three-dimensional  exponential  finite 
difference  equation  can  be  shown  to  be: 


IJ.k 


U.k 


Using  the  sub-time  Interval  concept,  Eq.  (91) 


becomes : 


Tn+1  xn  _ fi 
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♦ 1  Z-#  u. 
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where  m  Is  the  number  of  subtime  Intervals,  Q  Is  the 

dimensionless  time  step,  and  k  Is  the  dimensionless  drive  number 

given  by. 
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'U.k 
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Equation  (92)  will  be  used  for  all  Interior  nodes  In  Fig.  13.  This 
equation,  as  well  as  those  that  result  from  the  other  analysis,  will  be 
adjusted  along  the  Insulated  boundaries  to  take  Into  account  the 
boundary  condition  that  exists  there. 

The  next  method  to  be  applied  to  this  three-dimensional  case  will 
be  the  pure  explicit  method.  The  finite  difference  equation  for  this 
method  Is  given  by  [11  ]: 

IJ.k  =  T1  J,k  (  "  6Q)  +  Q  (TU1  J.k  +  T1-l  J.k  +  T1  ,J+l,k 


r n  Tn  *n  * 

1  J-l.k  +  'l  J.k+l  +  1  J,k-1  J 
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where  ft  *  ;  Ax  =  A y  =  az 

(AX)^ 

As  shown  In  Ref.  [11]  the  dimensionless  time  step  ft  must  be: 

n  <  l  (95) 

to  ensure  stability  of  the  method. 

The  last  numerical  technique  that  will  be  applied  Is  the  Method  of 
Douglas  [11].  This  method  Is  Implicit,  and  the  spatial  directions  are 
considered  sequentially  In  the  x,  y,  and  then  z  directions 
respectively.  The  Intermediate  temperatures  U  (found  from  x-dlrect 
sweep)  and  V  (found  from  y-dlrectlon  sweep)  are  used  to  calculate  the 
actual  temperature  field  variable  T  (found  from  z-dlrectlon  sweep). 

The  equations  that  are  solved  sequentially  are  presented  as  follows. 


U  -  T 

u1-1.k  1.1. k  1  .2 


(U1.1.k  +  T1J,k)  +  4y  (T1  J.k)  +  6z  (T1,j,k) 


Tn+1  Tn 
1.1. k  "  1.1. k 


*  4z  (T?.3,k)  <97> 

Tn+1  Tn  . 

1.1. k  J_  1,1. k  =  1  ^  J  k  J  1  ^ 

+  2  4z  (T1J,k  +  T1,1,k)  (98) 

where  the  finite  difference  operator  In  the  x-dlrectlon  ,  for  example, 
would  be: 

42  .  Hi.I.k  'i-l.l. k _ tllAA  (99) 

x  (ax)2 

Equations  (96),  (97),  and  (98)  must  be  solved  successively  because  the 
variable  U  Is  used  In  equation  (97)  to  find  V  and  so  on.  Since  the 
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method  operates  on  one  spatial  direction  at  a  time,  the  Thomas 
Algorithm  can  be  utilized.  In  the  case  of  finding  the  U  variable, 
the  y  and  z  nodal  positions  are  held  constant  for  all  the 

x-dlrectlon  nodal  positions  (Fig.  14).  This  process  Is  repeated  until 
all  y  and  z  nodal  values  for  the  x-dlrectlon  variable  U  are 

calculated.  This  procedure  Is  then  repeated  In  a  similar  way  for  the 
V  variable  and  then  finally  for  the  actual  temperature  field  variable. 

The  results  from  the  three  different,  three-dimensional  solution 
methods  are  shown  In  Table  VII.  The  exponential  finite  difference 
method  described  above  outperformed  the  pure  explicit  and  the  method  of 
Douglas  for  all  positions  as  shown  In  Table  VII. 

In  Ref.  [11]  nine  different  methods  to  solve  the  diffusion 
equation  In  th^ee  dimensions  were  Investigated.  The  method  of  Douglas 
was  the  preferred  method  because  of  Its  accurate  results  and  low 
computer  CPU  time.  In  that  study  the  pure  explicit  method  required  the 
lowest  amount  of  CPU  time  with  the  method  of  Douglas  requiring 
approximately  four  times  as  much.  In  the  present  study  all  three 
methods  were  run  on  two  different  mainframe  computers  to  Investigate 
how  well  these  three  methods  compared  In  CPU  times.  The  results  are 
shown  In  Table  VIII.  All  three  methods  were  exercised  for  the  same 
number  of  time  steps.  As  Indicated,  the  exponential  method  was 
approximately  three  times  faster  than  the  method  of  Douglas  but  still 
slower  than  pure  explicit  method.  From  these  results  It  could  be 
concluded  that  the  exponential  method  would  have  been  chosen  as  the 
preferred  method  had  It  been  used  In  competition  with  the  nine 
numerical  methods  as  described  In  Ref.  [11]. 


IV.  -  EXTENSION  OF  THE  EXPONENTIAL  FINITE  DIFFERENCE  METHOD 


TO  NONLINEAR  PARTIAL  DIFFERENTIAL  EQUATIONS 
The  exponential  finite  difference  method  will  now  be  applied  to 
two  different  nonlinear  problems.  The  problems  to  be  addressed  will  be 
the  viscous  Burger's  equation  and  the  boundary  layer  equations  (steady 
state  flow  over  a  flat  plate). 

The  viscous  Burger's  equation  Is  given  In  Ref.  [12]  as: 


a2U 


au  .  ..  au 

«  '  “  "  "  ax2 


(100) 


To  allow  application  of  the  exponential  method  to  this  equation,  the 
equation  must  be  first  linearized.  So  letting  U  =  A  =  constant,  for 
the  nonlinear,  term  and  rearranging  the  equation,  gives: 


3U 

at 


,  au  a2u 


(101) 


Assuming  a  product  solution  of  the  form 

U(x,t)  =  $(x)e(t) 

and  taking  the  appropriate  derivatives,  Eq.  (101)  becomes 


»  .  Ae  M  .  ve 


Division  by  gives: 


1  36  A  3<t>  v  3  d>  ,  , 

5  at  ■  -  ;  ax  *  ;  ~2  -  -  «  ■  constant 

o  X 


(102) 


It  can  be  seen  that  the  terms  are  now  separated.  As  has  been  shown 
earlier,  the  left-hand  side  of  Eq.  (102)  can  be  written  as: 
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U( x, t)  =  U(x,0)  exp  {-  ict }  (103) 

Now  returning  to  Eq.  (102)  and  examining  the  x-dependence,  we  have 

.2. 


A  3$  +  £  a  $ 
'  ♦  ax  4>  ax2 

Multiplying  both  sides  by  gives: 
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This  can  be  written  In  finite  difference  form  as 
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This  Is  used  to  replace  the  exponent  In  Eq.  (103),  thus 
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Equation  (105)  Is  the  exponential  finite  difference  equation  for  the 
viscous  Burger's  equation.  An  example  will  now  be  used  to  demonstrate 
the  method. 

An  exact  steady  state  solution  to  Burger's  equation  Is  available 
for  the  following  conditions 

U(0,t)  =  UQ 
U(L,t)  =  0 

The  steady-state  solution  was  given  as  [12]: 


U(x)  =  UQU1 


1  -  exp  (0  R 


_ -  j  t 

’  *  «p  (»  Re  (l  -  '))  J 


Rk  «  m 


where 


shown  to  be: 


The  results  obtained  by  applying  Eq.  (109)  and  the  conditions  In 

Eq.  (108)  are  compared  to  the  steady  state  exact  results  of  Eq.  (106) 

and  are  shown  In  Fig.  (15).  The  results  from  the  exponential  method 

were  nearly  the  same  as  the  exact  method.  The  exponential  method  was 

allowed  to  march  In  time  for  quite  a  number  of  steps  without  special 

2 

treatment  of  the  dimensionless  group  Atv/(Ax)  which  could  have 
been  altered  to  allow  convergence  to  the  actual  solution  In  less  time 
steps. 

Another  application  of  Burger's  equation  was  made  to  Investigate 
the  effect  of  the  diffusion  term.  The  results  for  the  variation  of  « 
over  four  orders  of  magnitude  are  shown  In  Fig.  (16)  for  the  same 
Instant  In  time.  At  the  two  lower  v  values,  the  total  range  of  the 
field  variable  takes  place  over  a  small  number  of  nodal  positions.  A 
better  approximation  could  be  made  for  these  cases  by  using  a  finer 
grid.  Also  Included  on  Fig.  (16)  Is  the  solution  of  Burger's  equation 
by  a  pure  explicit  technique.  For  the  value  of  v  chosen,  the 
solution  oscillates  around  the  predicted  solution  found  from  the 
exponential  method.  The  pure  explicit  solution  was  found  using  the 
same  number  of  nodes  and  the  same  dimensionless  step  size.  When  the 
solution  oscillates,  as  the  pure  explicit  solution  did,  the  resulting 
velocity  field  can  contain  physically  Impossible  values. 
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The  last  application  to  be  Investigated  will  be  for  the 
development  of  a  laminar  boundary  layer  on  a  flat  plate  (Fig.  (17)). 
In  Ref.  [10]  the  steady  state  formulation  Is  given  In  terms  of  the 
following  three  partial  differential  equations: 
for  continuity: 


ax  ay 


(HO) 


for  momentum: 


(J  iy.  +  y  1U 

ax  ay 


a2u 


(ill) 


ay 


for  energy: 


u  —  ♦  v  — 

ax  ay 


a2T 


(112) 


ay 


with  the  boundary  conditions: 

U(x,o)  =  0 


U(o,y)  -  U 


V(x,o)  =  0 
T( x,o)  =  0 


V(x,L)  »  0 
T(o,y)  =  T 


(113) 


v  and  a  are  the  momentum  and  thermal  dlf fusl vltles  respectively. 

Equations  (111)  and  (112)  can  be  solved  by  using  the  method 
presented  for  the  viscous  Burger's  equation.  The  only  difference  Is 
that  the  solution  will  march  In  the  x-dlrectlon  Instead  of  time. 
Keeping  this  procedure  In  mind,  results  of  the  separation  of  variables 
for  Eqs.  (Ill)  and  (112)  were  found  to  be: 
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The  continuity  equation  Is  written  as  £  1 2 J : 

u1+l  Ay  f.,1+1  ,,1  ..1+1  ,.1  \  .... 

-vj  -rsvuj  -  uj  *  Vi  -  Vi]  <m 

Equations  (114)  and  (115)  are  first  solved  using  a  spatial 
sub-increment  as  was  done  for  the  cases  when  time  was  the  marching 
direction  of  the  solution.  After  this,  the  continuity  Eq.  (116),  Is 
then  solved. 

The  results  of  this  application  are  shown  In  Fig.  (18)  for  a 
Prandtl  number  equal  to  0.72.  As  can  be  seen  the  thermal  boundary 
layer  was  outside  the  velocity  boundary  layer  as  would  be  expected. 

The  results  with  the  Prandtl  number  equal  to  0.72  were  compared  to  the 
exact  solution  as  presented  In  Ref.  [10].  A  downstream  position  was 
chosen  and  the  results  are  compared  In  Table  IX.  The  exponential 
method  results  were  In  good  agreement  with  the  exact  results. 


CONCLUDING  REMARKS 


In  conclusion,  an  exponential  finite  difference  technique  has  been 
extended  to  other  coordinates  systems  and  expanded  to  handle  problems 
In  two  and  three  dimensions.  The  method  has  direct  application  to 
linear  partial  differential  equations  such  as  the  diffusion  equation 
and  can  be  extended  to  solve  nonlinear  equations. 

The  method  was  applied  to  the  following  cases: 

(1)  One-dimensional,  unsteady  state  heat  transfer  In  cylindrical 
coordinates,  Infinite  and  finite  heat  transfer  coefficient. 

(2)  Two-  and  three-dimensional,  unsteady  state  heat  transfer  In 
Cartesian  coordinates. 

(3)  One  dimensional  heat  transfer,  with  temperature  varying  thermal 
conductivity. 

(4)  Developing  temperature  field  In  laminar  Couette  flow. 

(5)  Nonlinear  partial  differential  equations  (Burger's  equation 


and  boundary  layer  equations' 

The  exponential  finite  difference  method  predicted  the  field 
variable  with  a  higher  degree  of  accuracy  In  those  cases  examined  where 
the  exact  solution  was  available.  When  extended  to  three  dimensions, 
the. accuracy  was  still  higher  for  the  exponential  finite  difference  but 
the  computer  CPU  time  was  Increased.  When  the  exponential  method  was 
compared  to  other  numerical  techniques,  the  results  were  found  to  be 
very  comparlable. 
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In  conclusion,  the  results  predicted  for  the  exponential  finite 
difference  algorithm  for  the  cases  presented  In  this  study  demonstrated 
that: 

(1)  Field  variable  was  predicted  with  a  higher  degree  of  accuracy 
than  other  numerical  techniques  where  exact  solutions  exist. 

(2)  The  method  can  be  applied  to  linear  and  nonlinear  partial 
differential  equations  with  dependent  variables  that  can  be 
separated . 

(3)  The  stability  of  the  method  Is  the  same  as  that  of  pure 
explicit  methods,  where  the  sub-time  Interval  step  size 
determines  the  stability. 
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TABLE  I.  -  COMPARISON  Of  RESULTS  FOR  DIFFERENT  DIMENSIONLESS  TIME  STEP 
FOR  ONE- DIMENSIONAL  HEA1  1RANSFER  IN  CYLINDRICAL  COORDINATES  WITH 
INFINITE  HEAT  TRANSFER  COEFFICIENT  AT  THE  SURFACE.  INITIAL  AND 
BOUNDARY  CONDI! IONS  ARE: 

(h  -*  ®,  T(r.o)  *  T  .0,  T ( R , t)  -  0.0,  $1  -  a  At/(Ar)2,  a  =  1.0  M^/s. 

N  *  number  of  nodes,  m  *  number  of  sub-time  Intervals.) 

I  From  surface  I  N  =  TT  1  N  «  21  N  »  21  1  N  .  21  F  Ei 


t. 

From  surface 

N  =  11 

N  •  21 

sec 

r-dlstance 

m  =  4 

m  -  9 

(M) 

$1  =  1.0 

$1-1.0 

0.1 

0.1 

0.127004 

0.126/68 

.1 

1 .0 

.862431 

.852204 

.5 

.1 

.011959 

.011671 

.5 

1 .0 

.094334 

.090309 

.5 

Total 

Total 

50  steps 

200  steps 

-  — —  =  0  1 

M  ♦  1 

M  *  1  *  °'2 

0 . T  26819 
.853083 

.011680 

.090379 

Total 
TOO  steps 


.0  M?/s. 

21  1 

Exact 

9  j 

analysis 

5.0 

ref.  (9) 

_ 

0.126669 

.848368 

.011715 

.090652 

Total 
40  steps 


.011582 

.088895 


TABLE  II.  -  FINITE  HEAT  TRANSFER  COEFFICIENT 
CYLINDRICAL  COORDINATES  WITH  THE  FOLLOWING 
CONDI1 IONS:  T(r,o)  »  1.0, 

T  a  0,  FT  *  a  At/(  Ar) 


Exponential  finite 
difference  results 


T  Ime 

h/k 

R 

0.1 

i 

1 

0 

.2 

i 

1 

0 

.4 

i 

1 

0 

.1 

2 

1 

0 

.1 

5 

1 

0 

.9768  .9814 
.5702  .5976 
.8702  .8852 
.4132  .4441 
.6420  !  .6698 
.5009  !  .5285 
.9594  |  .9670 
.2558  i  .2777 
.9265  .9385 


N  =  1 1  N  «  21  N  .  21 

m  -  4  m«9  m  •  9 

$1*1.0  $1-5.0  $1-1.0 


.9797 

.5857 

.8780 

.4303 

.6563 

.5199 

.9643 


0.2669 

.9306 


.  r.  wVkVvWIG 
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TABLE  III.  -  COMPARISON  OF  EXPONENTIAL 
FINITE  DIFFERENCE  METHOO  IN  ONE- 
OIMENSIONAL  CYLINDRICAL  COORDINATES 
TO  THE  RESULTS  OF  REFERENCE  [4], 


[a  »  1.0,  At  -  1.0  sec,  a  At/Ar2  .  1.0. 
N  =  number  of  nodes,  m  =  number  of 
sub-intervals 


m 

Reference 

[4] 

sa 

N  .  19 
m  *  8 

ft  =  1  .0 

5 

D 

18 

0.77220 

0.773094 

0.772922 

1 

10 

.01449 

.011353 

.011951 

18 

.84661 

.846719 

.846811 

10 

.11595 

.112112 

.113523 

18 

.93546 

.935278 

.935521 

10 

.57722 

.575979 

.576198 

00 

18 

.99370 

.993686 

.993776 

10 

.95872 

.958596 

.959245 

TABLE  IV.  -  COMPARISON  OF  EXPONENTIAL  FINITE  DIFFERENCE 
METHOD  IN  TWO-OIMENSIONAL  CARTESIAN  COORDINATES 
TO  THE  ALTERNATING  OIRECTION  IMPLICIT  METHOD  [4]. 

[For  comparison  to  results  In  ref.  [4]  at  x  -  y  *  o.] 


Exact 

ADI 

$7  »  1.0,  N  -  11, 

ft 

«  1.0,  N  .  21, 

[4] 

At  *  0.01,  m  ■  4 

At 

=  0.0025,  m  .  9 

0.1 

0.09883 

0.09333 

0.09829 

0.09924 

.2 

.40354 

.40354 

.40256 

.40354 

.3 

.63179 

.63224 

.63080 

.63166 

.4 

.7  '486 

.77532 

.77403 

.77472 

.5 

.86252 

.86283 

.86187 

.86240 

.6 

.9160/ 

.91624 

.91559 

.91597 

.7 

.9487/ 

.94886 

.94841 

.94869 

.7 

(Total  of 

(Total  of 

70  steps) 

280  steps) 

n  n  - 

_  _ 

m  ♦  1  =  °'2 

-  i  *0.1 

m  +  1 
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TABLE  V.  -  STEADY  STATE  HEAT  TRANSFER  IN  TWO-DIMENSIONS 

Comparison  of  exponential  finite  difference  technique  to  a  Gauss-Seldel 
technique  for  the  solution  of  Laplace's  equation. 


9  by  9  Grid 


Gauss-Seldel 
88  Iteratlcns(a) 


0.000 

1.000 

0.0 

0.0 

.125 

.875 

1.7413 

1.7414 

.250 

.750 

6.8946 

6.8949 

.375 

.625 

15.0330 

15.0335 

.500 

.500 

24.9999 

25.0004 

.625 

.375 

34.9667 

34.9672 

.750 

.250 

43.1052 

43.1055 

.875 

.125 

,  48.2587 

48.2588 

1.000 

0.000 

100.000 

100.000 

aFrom  reference  [4]. 


5  by  5  Grid 


Exponential 
finite 
difference 
*0  Iterations 


Gauss-Seldel 
22  Iteratlons(a) 


0.0 

7.1428 

25.0000 

42.8571 

100.000 


Exponential 
finite 
difference 
20  Iterations 


7.1430 

25.0003 

42.8573 


100.000 


TABLE  VI.  -  COMPARISON  OF  EXPONENTIAL,  PURE-EXPLICIT,  AND 
IMPLICIT  FINITE  DIFFERENCE  METHODS  FOR  ONE-DIMENSIONAL, 
UNSTEADY-STATE  HEAT  TRANSFER  WITH  TEMPERATURE  VARYING 
THERMAL  CONDUCTIVITY  AT  THE  CENTER  OF  THE  SLAB 
(X(T)  .  1.0  ♦  8(T);  8  -  0.01.1 


L*  •  *V 


wy ' 


i  Lk  ,  •  *  •  at  »  1  $,♦  . .  ■  ,  -  .  »  *,«  1,1  Mj,*  »,*  ♦  3,?  »*•  i.  '  Cl  -»>,•  »>  » .’  <;'  e.»  4.  »t  .4,'  U  .««  _u  V^i^j  •!  «  1  ,~t  »  |  *  t  r*f  «*♦  |’I_| 
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TABLE  VII.  -  COMPARISON  OF  THREE  DIFFERENT,  THREE-DIMENSIONAL  UNSTEADY  STATE 

HEAT  TRANSFER  SOLUTIONS 


[T(x,y,z,o)  -  1.0;  T(x,y,L,t)  -  T(x,L,z,t)  .  T(L.y.z.t)  »  0;  aT/ax(o,y,z,t)  • 
dT/ay(x,o,z,t)  *  aT/az(x,y,o,t)  «  0;  N  ■  number  of  nodes  In  x,  y,  and  z 
directions;  Q  «  a  At/(Ax)2  and  Ax  ■  ay  ■  Az.] 


Elapsed 

time, 

sec 

Position  from 
center  along 
diagonal 
x  •  y  •  z 

Exact 

analysis 

result, 

°C 

Exponential  finite 
difference  results, 
°C 

N  >  11 ,  m  •  4, 
a  .  0.75 

(a) 

Pure  explicit 
finite  difference 
results, 

•c 

Q  •  0.15,  N  .  11 

(a) 

Method  of  Douglas 
Douglas  finite 
difference  results, 

•c 

Q  •  0.15.  N  •  11 

(a) 

0.09 

0.0 

0.693490 

0.892237  (0.14) 

0.889437  (0.45) 

0.886760  (0.75) 

.5 

.440712 

.440650  (  .014) 

.435058  (1.28) 

.439665  (  .24) 

.9 

.006491 

.006484  {  .11) 

.006319  (2.65) 

.006510  (-.29) 

.15 

.645469 

.645209  (  .04) 

.640025  (  .84) 

.641484  (  .62) 

.5 

.253065 

.253286  (-.09) 

.250102  (1.17) 

.252641  (  .17) 

.9 

.003015 

.003022  (-.23) 

.002970  (1.49) 

.003023  (-.27) 

aAccuracy  percent. 


TABLE  VIII.  -  COMPARISON  OF  C.P.U.  TIME  ON  TWO 
DIFFERENT  MAINFRAMES  FOR  THREE  DIFFERENT 
THREE-DIMENSIONAL  FINITE  DIFFERENCE  METHODS 


[One-hundred  time  steps  for  each  method.) 


Computer 

Exponential* 

Method  of 

Pure-explicit 

method. 

Douglas, 

method, 

sec 

sec 

sec 

CRAY-XMP 

0.2778 

0.955 

0.0627 

IBM-3033 

5.4 

12.6 

1 .0 

aBased  on  the  number  of  sub-time  Intervals 
equal  to  100. 


TABLE  IX.  -  COMPARISON  OF  EXPONENTIAL 
FINITE  DIFFERENCE  METHOD  TO  EXACT 
RESULTS  OF  BOUNOARY  LAYER  EQUATION 
[10]  FOR  THE  VELOCITY  PROFILE  AT 
ONE  OOWNSTREAM  LOCATION. 


[Distance  downstream  x  -  500  cm, 
*  »  0.0072  cm2/s . ] 


Distance 
perpendicular 
to  plate, 
y  (cm) 

Exact 

result 

[10] 

Exponential 
method  result, 

N  *  21  m  *  8 

1.0 

0.17 

0.17428 

.34 

.34643 

.51 

.51020 

K] 

.65658 

R9 

.77684 

.87 

.86636 

.93 

.92638 

8.0 

.96 

.96265 
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i.j  -  SPACIAL  GRID 
n  -  Tilt  GRID 
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1st  2nd 

TIME  SUB-  TINE  SUB- 

I MTERVA1.  INTERVAL 


n 

NODE  i-1  • 


FINE  STEP 

n»1/3  n*2/3  n+1 

’  T 

Ax 


101  Al  DIFNNMONUSS  MIN  SUP 

FIGURl  7.  COflPUTA!  IONA!  GRID  fOR  7  Hrt  SDR  INIIRVAIS  CM  -  7).  CARI 
TLS1AN  COORDINAII  S.  (SHOWN  FOR  1-DIXNSION). 


PROGRAM 

.VARIABLES 


l=i.. .mot  y 
i  £  ntot 

I  1...N  \ 

i  1...N  / 


Pl  -I  °  o 
v,M=v..i 


K*1.2. .NS1 


J=1.2...N 

M.2...N 


I  m  .VTi*M4VTM.l*VT».l*i4¥Ti.i-rWTi. 


- */  "1-2 . ■ 


I— T  VT4  j-VTj  t  <£XPaM((/MS1> 


■M.7....N 


' — i  p'-rp»-.,"i-i 


It. 7 . N 

J ■  1.7 . N 


Vl  -j 


FIGURE  J.  -  FLOW  DIAGRAM  FOR  2-DIMENSIONAL  EXPONENT IAE 
FINITE  DIFFERENCE  ALGORITHM. 

NTOT  -  TOTAL  NUMBER  OF  TIME  STEPS 
P1>(  -  SUM  OF  DIMENSIONLtSS  DRIVE  NUMBERS 
VT j  :  -  field  variable  during  time  siep 

M  ’  -  DIMENSION!  ESS  DRIVE  NUMBER 
NS1  -  NUMBER  OF  TIME  SUB- I NTERVALS 
001 

0  -  DIMENSIONLESS  TIME 


5»Ss 


UNSTABLE 

OSCILLATIONS 

FIGURE  A.  -  EFFECT  OF  NONDINENSIONAL  TINE  STEP  SIZE  ON  1  NODE  HOOEl 
SOLUTION.  (FROM  REF.  (21). 


THERMAL  CONDUCTIVITY 


.s 

CROSS-STRf AH  POSITION,  (m> 


1.0 


HF-URl  17.  WVHOPINC  TEHPt  RAIURE  FlflD  IN  IABINAR  COUETTE 
1 1  ON.  SHOWING  THE  EPFTCT  OF  F1UID  VTIOCITY  'O’:  (l  1.0  M>. 


v*sVI\v 
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I  IGURl  16.  iXPOWNIAt  i  INI1I  DllftRINCL  RLSULIS  FOR  VARY¬ 
ING  KIM  Ml  1C  VISCOSITY.  ALL  VHOCIIItS  ARC  SHOWN  FOR 


1  1.0  StCOHDS  WITH;  N  ?1.  H  =  8.  —  i  --  4.0. 

(AKU 


U(o.t)  *  1.0 
IIU.1)  --  0.0 


(AX) 


I  I  MINI  I/. 


01  lit  K  CONDI  HONS:  U(x,o)  0.  V(x.O)  0 
V(x.l  )  0 

WXINDAKY  I  AYt  K  IXVilomAI  Al  ONI,  A  COM  FI)  II  A I  IIAll. 


100 


300 


200 


400  WO  0  100  200  300  4005 

X  -  DISTANCE  DOWN  THE  FIAT  PLATE,  cm 


riGURt  18.  -  EXPONENTIAL  FINIIE  DIFFERENCE  RESUE LS  FOR  BOUNDARY  LAYtR  tOUATIONS.  WITH  CONDITIONS  U(o.y)  *  1.0. 
IKx.o)  0.  V(X.O)  '  0.  V(x.l)  0.  T(x.o)  0.0.  Ho.y)  1  1.0;  V  -  0.00/2  CM?/SIC;  0  -  0.01  CM?/SEC. 


m 


APPENDIX 

This  appendix  contains  all  the  computer  programs  mentioned  In  this 
report.  A  computer  program  variable  list  Is  also  contained  with  a 
description  of  their  use,  and  a  program  number  to  refer  to  the  programs 
that  they  are  contained  In. 

Each  of  these  programs  was  written  to  be  run  In  an  Interactive 
mode  with  the  mainframe  computer.  The  only  cases  run  differently  were 
for  the  three-dimensional  unsteady  state  heat  transfer  cases  that  were 
run  In  batch  mode  on  the  Cray  X-MP. 

The  program  structure  Is  as  follows.  A  main  program  Is  used  to 
describe  the  necessary  parameters  and  for  asserting  the  proper  boundary 
conditions.  The  main  program  then  calls  the  subroutine  where  the 
actual  finite  difference  methods  are  exercised  and  the  results  are  then 
printed. 
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COMPUTER  PROGRAM  LIST 


Program 

number 

Program 

name 

Program  function 

1 

SOURCE. EFOCYL 

One-dimension,  unsteady  state,  cylindrical 
coordinates.  Infinite  and  finite  heat 
transfer  coefficient 

2 

SOURCE. EFD2D 

Two-dimensional  Cartesian  coordinates, 
unsteady  state  heat  transfer 

3 

SOURCE. LAPLAC 

Two-dimensional  Laplace's  equation 

4 

SOURCE. EFOVAR 

One- dimensional  unsteady  state  heat 

conduction,  varying  thermal  conductivity, 
exponential  finite  difference  method 

5 

SOURCE. EXPVAR 

One-dimensional,  unsteady  state  heat 
conduction,  varying  thermal  conductivity, 
explicit  finite  difference  method 

6 

SOURCE. IMPVAR 

One-dimensional,  unsteady  state  heat 
transfer,  varying  thermal  conductivity, 
Implicit  finite  difference  method 

7 

SOURCE. COUE 

One-dimensional,  developing  temperature  field 
In  laminar  couette  flow 

8 

SOURCE. EX30 

Exact  analysis,  three-dimensional  heat 
transfer  In  a  cube 

9 

SOURCE. EF03D 

Three-dimensional  unsteady  state  heat 
transfer  In  a  cube  using  exponential  finite 
difference  method 

10 

SOURCE. EXPL3D 

Three-dimensional  unsteady  state  heat 
transfer  In  a  cube  using  explicit  finite 
difference  method 

11 

SOURCE. OOUGLA 

Three-dimensional  unsteady  state  heat 
transfer  In  a  cube  using  the  method  of 
Douglas 

12 

SOURCE. BURGER 

Exponential  solution  of  nonlinear  viscous 
Burger's  equation 

13 

SOURCE. EXBURG 

Pure  explicit  solution  of  nonlinear  viscous 
Burger's  equation 

14 


SOURCE .NONBOU  Exponential  Method  of  solution  for  boundary 
layer  equations  for  flow  over  a  flat  plate 
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COMPUTER  PROGRAM  VARIABLE  LIST 


Program 

variable 

name 

Programs 
used  In 

Variable  Description 

N 

1-14 

Number  of  nodes 

NS 

1-12,14 

Number  of  time  sub  Intervals 

NTOT 

1-14 

Total  number  of  time  or  spaclal  steps 

TSI 

1-13 

Dimensionless  time  step  Increment 

T 

1-7,9,10,12,13 

Total  elapsed  time  or  spatial  distance 
between  steps 

DL 

1 

Radial  distance  between  adjacent  nodes 

R 

1 

Radial  length 

I  PR 

1-14 

Number  of  steps  between  output  of  results 

NB 

1 

Heat  transfer  boundary  condition  flag 

V 

1-14 

Dependent  variable 

HK 

1 

Convection  heat  transfer  coefficient 
divided  by  thermal  conductivity 

VM 

1 

Dependent  variable  value  outside  of 
solid  In  the  surrounding  medium 

M 

1-10,12 

Dimensionless  drive  number 

P 

1-10,12 

Sum  of  the  drive  numbers 

VT 

1-14 

Dependent  variable  value  during 
sub-time  Interval 

B 

1 

Blot  number 

THE 

1-3,12-14 

Variable  used  for  the  output  of  results 

TIME 

1-6,8-13 

Total  elapsed  time  of  the  solution  at 
the  current  output 

UMAX 

1-14 

Output  counter 

ET1 

3 

Accuracy  desired  In  solution  of 

Laplace's  equation 

Program 

variable 

name 

Program 
numbers 
used  In 

Description 

DELV 

3 

Difference  In  dependent  variable  value  from 
one  time  step  to  the  next 

ETA 

3 

Sum  of  the  absolute  value  of  the  differences 
found  In  DELV 

KR 

4-6 

Reference  thermal  conductivity 

BETA 

4-6 

Slope  of  thermal  conductivity  variation  with 
temperature 

K0,K 

4 

Thermal  conductivity  at  the  total  time  step 
Interval  or  sub-time  Interval  respectively 

THE 

4 

Klrchoff  transformation  variable  (used  In 
exponential  finite  difference  program  with 
varying  thermal  conductivity) 

KAPPA 

4 

All  values  known  from  the  last  time  step 
Increment  and  used  to  solve  the  quadratic 
equation  that  results  In  the  exponential 
finite  difference  solution  with  temperature 
varying  thermal  conductivity 

DERSQR 

5 

Absolute  value  of  velocity  difference  found  1i 
evaluation  of  velocity  gradient 

OHEGA 

6 

Same  as  nondlmenslonal  time  step 

A,B,C,D 

6,11 

Coefficient  used  In  trldlagonal  matrix 
algorithm. 

KAP.GAM 

6 

Variables  used  to  determine  A,B,C 

JETA, GAMMA 

6 

Variables  used  In  Thomas,  trldlagonal 
algorithm 

TS 

6 

Same  as  total  elapsed  time 

T 

6 

Dependent  variable 

V 

6 

Solution  vector  tri-diagonal  algorithm 

SP 

7 

Maximum  width  divided  by  maximum  velocity 

FL 

7 

Parameter  based  on  position  In  flow 

OIST 

7 

Serves  same  function  as  time  for  unsteady 
state  problem 

n 

Program 

variable 

name 

Program 
numbers 
used  In 

Description 

PI , PI2, PI3 

8 

it,  ir^,  *3 

NOOES 

8 

Same  an  N 

TO 

8 

Initial  temperature 

TI 

8 

Surface  temperature,  t  >  0 

ALFA 

8 

Thermal  dlffusivlty 

TNEW 

8 

Exact  temperature  at  a  x,  y,  and  z 
location  after  elapsed  time  t  has  occurred 

T 

11 

Dependent  variable 

OELX.DELY, DELZ 

11 

Part  of  the  central  difference  operator 

u.v 

11 

Variables  used  to  sweep  preliminary  solution 

In  the  x  then  y  directions  respectively 

H 

11 

Used  as  an  array  to  contain  the  known 
quantities  used  In  the  trldlagonal  algorithm 

RNU 

12-14 

Kinematic  viscosity 

DX 

14 

Step  length  In  flat  plate  direction 

RAL 

14 

Thermal  dlffusivlty 

YMAX 

14 

Maximum  distance  perpendicular  to  flat  plate 

DY 

14 

Step  length  perpendicular  to  the  flat  plate 

u.v 

14 

Dependent  variables  (velocity)  In  x  and  y 
directions  respectively 

T 

14 

Temperature  field  variable 

MU,  Ml 

14 

Drive  numbers  for  velocity  In  x-dlrectlon 
and  temperature  field 

PU.PT 

14 

Sum  of  drive  numbers  for  MU,  MT 

Ul.TT.VT 

14 

x-dlrectlon  velocity,  temperature,  and 
y-dlrectlon  velocity  on  subinterval 

UT1 

14 

Temporary  U-dlrectlon  velocity  field 

IS.TS1 

14 

Dimensionless  time  step  for  temperature  and 
x-dlrectlon  velocity  respectively 

nnn  noo  nnnnnnnnnnnonoi 


'1i  4'^  fa  4  'a  ■  4  **  >.»  >  Q  <■»  »> 
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WRITTEN  BY  R.F.  HANDSCHUH 

SOURCE. EFDCYL 

*****  PROGRAM  #1  ***** 


THIS  PROGRAM  IS  TO  BE  USED  AS  THE  STARTING  POINT  FOR  INVESTIGATING 
THE  EXPONENTIAL  FINITE  DIFFERENCE  ALGORYTHM.  THIS  METHOD  WAS  INTRODUCED 
BY  M.C.  BHATTACHARYA.  THIS  SOURCE  CODE  IS  FOR  CYLINDRICAL 
COORDINATES,  UNSTEADY-STATE  HEAT  CONDUCTION,  1  DIMENSION. 


IMPLICIT  REAL#8! A-H.O-Z) 

REAL*8  V(100) ,R(100) 

INPUT  THE  PROGRAM  DATA 

WR1TE(6 , 15) 

5  FORMAT! IX, ’NUMBER  OF  NODES=N  13’/) 

READ( 9, 10 )N 
10  FORMAT! 13) 

WRITE! 6, 12) 

12  FORMAT! IX, ’NUMBER  OF  TIME  SUB  INTERVALS1  NS  13’) 

READ! 6 , 1 3 )NS 

13  •  FORMAT! 13) 

WRITE!6 , 16 ) 

16  FORMAT! IX, ’TOTAL  NUMBER  OF  TIME  STEPS*  NTOT  13’) 

READ! 9 ,21 )NTOT 

21  FORMAT! 13) 

WRITE! 6, 24) 

24  FORMAT! IX, ’INPUT  TIME*THERMAL  DIFFUSIVITY  /  RAD  INT  SQUARED  F5.3’) 
READ! 9,25 )TSI 

25  FORMAT !F5. 3) 

WRITE! 6, 22) 

22  FORMAT! IX, ’TOTAL  TIME  OF  ONE  TIME  STEP*  T  F7.5’) 

READ! 9 , 23 )T 

23  FORMAT !F7. 4) 

WRITE! 6, 26) 

26  FORMAT!  IX, ’INPUT  RADIAL  INTERVAL  I,ENGTH=DL  F5.3’) 

READ! 9,27)DL 

27  FORMAT! F5. 3) 

R !  1 )  —  1 . 

DO  28  1=2, N 
IM1=I-1 

28  R!I)=R(IM1)-DL 
WRITE! 6, 32) 

32  FORMAT! IX, ’INPUT  NUMBER  OF  TIME  STEPS  BEFORE  PRINTING  13’) 

READ! 9,33 )IPR 

33  FORMAT! 13 ) 

DETERMINE  THE  TYPE  OF  BOUNDRY  CONDITION,  THEN  SET  VALUES 


noon  non  ono  oon  ooo  non 
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WRITE ( 6,14) 

14  FORMAT ( IX, 'INPUT  HEAT  TRANSFER  B.C.  0  -  INFINITE  1  -  FINITE'/) 
READ( 9 , 17 )NB 

17  FORMAT (II ) 

IF(NB.EQ.l)  GO  TO  100 
V(1)  =  0  . 

DO  30  1=2, N 

30  V(I)=1.0 

CALL  EXP  FIN  DIF  FOR  INFINITE  HEAT  TRANSFER  COEFFICENT 

CALL  EFDIHC( N , NS , NTOT , TSI , V , T , R , DL , IPR ) 

GO  TO  101 

100  CONTINUE 
WRITE( 6,31) 

31  FORMAT ( IX, 'INPUT  HEAT  COEF  /  TERM  COND  F5 . 3 ' ) 

READC  9 , 34 )HK 

34  FORMAT (F5. 3) 

DO  40  1=1, N 
40  V(I)=1. 

VM=0 . 

CALL  EXP  FIN  DIF  FOR  FINITE  HEAT  TRANSFER  COEFFICENT 

CALL  EFDFHC ( HK . N , VM , DL , NS , V , NTOT , TSI , T , R , IPR ) 

101  CONTINUE 
STOP 
END 

SUBROUTINE  EFDIHC 

SUBROUTINE  EFDIHC ( N . NS , NTOT , TSI , V , T , R , DL , IPR ) 

FOR  INFINITE  HEAT  TRANSFER  COEFFICENT 
IMPLICIT  REAL#8(A-H,0-Z) 

REAL*8  VT(10Q),V(100),M(100),P(100),R(100), THEC 100) 
TS=TSI/DFLOAT(NS+l ) 

WRITE ( 6,21)(R(I),I=1,N) 

21  FORMAT ( 1X,11(F6.3,2X)) 

WRITE ( 6 , 22 )DL , TSI 

22  FORMAT ( 1X,2(F6.3,2X)) 

N1=N-1 

NS1=NS+1 

BEGIN  MAIN  TIME  STEP  LOOP 
DO  20  L=1 , NTOT 

ZERO  DRIVE  NUMBERS  AND  SET  TEMPORARY  VARIABLES  EQUAL  TO  THE 
LAST  TOTAL  TIME  STEP  VALUES 

DO  15  1=1, N 

15  P(I)=0. 
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DO  10  1=1, H 
10  VT(I)=V(I) 

SUB  TIME  INTERVAL 

DO  30  K=1 , NS1 

CALCULATE  THE  DRIVE  NUMBERS 

DO  40  1=2, N1 

IM1=I-1 

IP1=I+1 

MCI)=(2.«VT(I)-VT(IM1)-VT(IM ) )/VT(I) 

40  MCI)=MCI)+DL#C VT(IP1 )-VTCIMl ) 1/C VTC I )*2 . «RCI) ) 
MCN)=2.*(VT(N)-VT(N1))/VT(N) 

CALCULATE  THE  DEPENDENT  VARIABLE  ON  THE  SUB-INTERVAL  LEVEL 
DO  50  11=2, N 

50  VT(I1)=VT(I1 )*DEXP( -TS#M( II ) ) 

SUM  THE  DRIVE  NUMBERS 

DO  60  1=2, N 
60  PCI)=PCI)+M(I) 

30  CONTINUE 

CALCULATE  THE  DEPENDENT  VARIABLE  ON  THE  NEXT  COMPLETE  TIME  STEP 
DO  70  1=1, N 

V(I)=VCI)#DEXPC-TS*PCI) ) 

70  CONTINUE 
ITMAX=ITMAX+1 

PRINT  THE  RESULTS 

IFC ITMAX . LT . IPRJGO  TO  20 
DO  71  1=1, N 

71  THECI)=VCI) 

ITMAX=0 
WRITEC  6,5) 

5  FORMAT ( / ) 

WRITEC  6 , 31 )L 

31  FORMAT ( IX, 'TIME  STEP  NUMBER=',I3) 

TIME=T*DFLOATCL) 

WRITEC  6 , 32 )TIME 

32  FORMATC1X, 'ELAPSED  TIME= ' ,F10. 4, '  SECONDS') 

IFCN.GT. ll)N21=N/2 

IFCN.GT. 11 )GO  TO  81 
WRITEC6 ,82) CTHECI) ,1=1 ,N) 

82  FORMAT ( 11C2X.F8.6)) 

GO  TO  84 
81  CONTINUE 

WRITEC 6 ,82) CTHECI) ,1=1 ,N21 ) 
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NNS-N21+1 

WRITE ( 6,82) (THE(I) , I=NNS,N) 

84  CONTINUE 
20  CONTINUE 
RETURN 
END 

SUBROUTINE  EFDFHC 

SUBROUTINE  EFDFHC  C  HK , N , VM , DL , NS , V , NTOT , TSI , T , R , IPR ) 

IMPLICIT  REAL*8  (A-H.O-Z) 

REAL*8  VT(IOO) ,V(100) ,M(100).P(100),R(100) ,THEC100) 

B=HK#DL 

TS=TSI/DFLOAT( NS+1 ) 

NS1 =NS+ 1 

N1=N-1 

N2=N-2 

BEGIN  THE  TOTAL  TIME  STEP  LOOP 
DO  20  L=1 ,  NTOT 

ZERO  THE  DRIVE  NUMBERS  AND  SET  THE  TEMPOARY  DEPENDET  VARIABLES 
EQUAL  TO  THE  LAST  COMPLETE  STEP  VALUES 

DO  15  1=1, N 
15  P(I)=0. 

DO  10  1=1, N 
10  VT(IJ=V(I) 

SUB  TIME  INTERVAL 

DO  30  K=1 , NS1 

CALCULATE  THE  DRIVE  NUMBERS 

M(l)=-(2.*VT(2)-(2.+2.*B)*VT(l)+2.«B*VM)/VT(l) 

M(1)=M(1)-DL*B*(VT(1)-VM)/<RC1)*VT(1)) 

DO  40  1=2, N1 
IM1 =1- 1 
IP1=I+1 

M(I)=(2.#VT(I)-VT(IM1)-VT(IP1 ) )/VT(I) 

40  MCI)=M(I)+DL*(VT(IP1 )-VT(IMl ) )/(VT(I)*2 ,#R(I) ) 
M(N)=2.*(VT(N)-VT(N1) )/VT(N) 

CALCULATE  THE  DEPENDENT  VARIABLE  ON  THE  SUB-INTERVAL  LEVEL 
DO  50  11=1, N 

50  VT(I1)=VT(I1 )#DEXP( -TS*M( II ) ) 

SUM  THE  DRIVE  NUMBERS 

DO  60  1=1, N 
60  P(I)=P(I)+M(I) 
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CONTINUE 


CALCULATE  THE  DEPENDENT  VARIABLE  ON  THE  NEXT  COMPLETE  TIME  STEP 


DO  70  1=1, N 

VCI)=VCI)*DEXPC-TS*P(I) ) 

CONTINUE 

ITMAX=ITMAX+ 1 

IF C ITMAX . LT . IPR ) GO  TO  20 


PRINT  THE  RESULTS 


DO  71  1=1, N 
THECI)=VCI) 

ITMAX=0 
WRITEC  6,5) 

FORMATC/) 

WRITEC 6, 31 )L 

FORMATC IX, 'TIME  STEP  NUMBERS,  13) 
TIME=T*DFLOATCL) 

WRITEC  6 , 32 )TIME 

FORMATC IX. ’ELAPSED  TIME= ' , FI 0 . 4 , ’  SECONDS') 
IFCN.GT. 11 )N21=N/2 
IFCN.GT.1DGO  TO  81 
WRITEC  6,82) CTHEC I) , 1=1 ,N ) 

FORMATC 11C2X,F8  &)) 

GO  TO  84 
CONTINUE 

WRITEC  6,82) CTHEC I) ,1=1 ,N21 ) 

NNS=N21+1 

WRITEC  6,82) CTHECI) ,I=NNS, N) 

CONTINUE 

CONTINUE 

RETURN 

END 
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SOURCE. EFD2D 

WRITTEN  BY  R.F.  HANDSCHUH 
*****  PROGRAM  #2  ***** 

THIS  PROGRAM  IS  FOR  2-DIMENSIONAL  CARTESIAN  COORDINATES 
UNSTEADY  STATE  HEAT  TRANSFER.  THE  METHOD  OF  SOLUTION  IS  THE 
EXPONENTIAL  FINITE  DIFFERENCE  ALGORYTHM.  THIS  PARTICULAR  PROGRAM  IS 
FOR  INFINITE  HEAT  TRANSFER  COEFFICENT  AT  THE  EXPOSED  SURFACES 
AT  X=Y=1 . 0  FOR  X  =  Y=0  THE  SURFACE  IS  CONSIDERED  TO  BE 
PERFECTLY  INSULATED. 


IMPLICIT  REAL*8(A-H,0-Z) 

REAL*8  V(25,25) 

INPUT  PROGRAM  DATA 

WRITE! 6, 15) 

15  FORMAT! IX, ‘NUMBER  OF  NODES=N  13’/) 

READ! 9 , 1 0 )N 

10  FORMAT (13) 

WRITEC  6,12) 

12  FORMAT! IX, ‘NUMBER  OF  TIME  SUB  INTERVALS*  NS  13') 

•  READ! 9, 13)NS 

13  FORMAT! 13) 

WRITE! 6, 16) 

16  FORMAT! IX, ‘TOTAL  NUMBER  OF  TIME  STEPS*  NTOT  I3‘) 

READ! 9,21 )NTOT 

21  FORMAT! 13 ) 

WRITE (6, 24) 

24  FORMAT! IX, 'INPUT  TIME*THERMAL  DIFFUSIVITY  /  LENGTH  SQUARED  F5 . 3 ' ) 
READ! 9 , 25 )TSI 

25  FORMAT(F5 . 3) 

WRITE(6 ,22) 

22  FORMAT! IX, 'TOTAL  TIME  OF  ONE  TIME  STEP*  T  F6.4') 

READ! 9 , 23 )T 

23  FORMAT! F6. 4) 

WRITE! 6 ,26 ) 

26  FORMAT! IX, ‘NUMBER  OF  TIME  STEPS  BEFORE  PRINTING*  I3‘) 

READ! 9 , 27 )IPR 

27  FORMAT! 13) 

WRITE! 6 , 250 )N, NS, NTOT 

250  F0RMAT(1X, *#  OF  NODES*' ,I3,2X, '*  OF  SUB-TIME-INT*' ,I3,2X, 

*'#  OF  TIME  STEPS* ',13) 

WRITE! 6 , 25 1 )TSI , T 

251  FORMAT! IX, ' (TIME*THER  DIFF)/LENGTH  SQUARED* ' F5 . 3 , 2X , 

* ' TIME  STEP  LENGTH*’ ,F6 .4/) 

N1  =  N-1 

INITIALIZE  THE  BOUNDRY  CONDITIONS 
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DO  30  1=1, N1 
DO  30  J=1,N1 
30  V(I,J)=1.0 

C 

I=N 

DO  50  J=1 ,N 

50  V(I,J)=0.0 
J=N 

DO  51  1=1, H 

51  V(I,J)=0.0 

CALL  EXP  FIN  DIF  FOR  INFINITE  HEAT  TRANSFER  COEFFICENT 

CALL  EFDIHC( N , NS , NTOT , TSI , V , T , IPR ) 

STOP 

END 

SUBROUTINE  EFDIHC 

SUBROUTINE  EFDIHC( N . NS , NTOT , TSI , V , T , IPR ) 

FOR  INFINITE  HEAT  TRANSFER  COEFFICENT 
IMPLICIT  REAL#8( A-H.O-Z) 

REAL#8  VT(25,25),V(2j,25),M(25,25),P(25,25), THE ( 25 , 25 ) 

PRINT  HEADING 
MRITEC  6 , 222 ) 

222  FORMAT(lX, ************  SOURCE. EFD2D  ##*#*#*##*#» •// ) 

TS=TSI/DFLOAT(NS+l) 

N1=N-1 

NS1=NSU 

BEGIN  MAIN  TIME  STEP  LOOP 
DO  20  L= 1 , NTOT 

ZERO  THE  DRIVE  NUMBERS  AND  SET  TEMPORARY  DEPENDENT  VARIABLE 
EQUAL  TO  THE  LAST  FULL  TIME  STEP  VALUE 

DO  15  J=1,N 
DO  15  1=1, N 
15  P(I,J)=0. 

DO  10  J= 1 , N 
DO  10  1=1, N 
10  VT(I,J)=V(I,J) 

SUB  TIME  INTERVAL 

DO  30  K=1 ,NS1 

CALCULATE  THE  DRIVE  NUMBERS 
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DO  41  J=2  ,  N1 
JM1 =J- 1 
JP1=J+1 
DO  40  1=2, HI 
1(11=1-1 
IP1=I+1 

M(I,J)=(VT(IP1,J)+VT(IM1,J)+VT(I,JP1)+VT(I,JM1)-4.*VT(I,J)) 

M(I, J)/V7(I, J) 

CONTINUE 

INSULATED  BOUNDRY  ALONG  X-AXIS 
J=1 

JP1=J+1 
DO  42  1=2, N1 
IP1=I+1 
IM1=I-1 

MCI,J)=(VT(IP1 , J) +VTCIM1 , J)+2 . #VT(I, JP1 )-4 . #VT(I, J) )/VT(I, J) 
INSULATED  BOUNDRY  ALONG  Y-AXIS 
1=1 

IP1=I+1 
DO  43  J=2,N1 
JP1=J+1 
JM1=J-1 

M(I,J)=(2. #VT(IP1 , J)+VT(I, JP1 )+VT(I, JM1 )-4 . *VT(I» J) )/VT(I , J) 
CORNER  AT  ORIGIN 

M(l,l)=(2. *VT( 1 , 2 ) +2 . *VT(  2 , 1 )-4 . #VT(  1 , 1 ) )/VT( 1,1) 

CALCULATE  THE  DEPENDENT  VARIABLE  ON  THE  SUB-INTERVAL  LEVEL 

DO  50  11=1, N1 
DO  50  J1 =1 , N1 

VT(I1,J1)=VT(I1,J1)#DEXP(TS*M(I1,J1)) 

SUM  THE  DRIVE  NUMBERS 

DO  60  1=1, N1 
DO  60  J=1,N1 
P(I.J)=P(I, J)+M(I,J) 

CONTINUE 

CALCULATE  THE  DEPENDENT  VARIABLE  AT  THE  NEXT  COMPLETE  TIME  STEP 


DO  70  J=1,N 
DO  70  1=1, N 

V(I, J)=V(I, J)#DEXP(TS*P(I,J)) 

CONTINUE 

ITMAX-ITMAX+1 

IFC  UMAX  .  LT .  IPR  )GO  TO  20 


HI 


C  PRINT  THE  RESULTS 

C 

WRITE( 6 , 5 ) 

5  FORMATC/) 

WRITE( 6 , 31  )L 

31  FORMATC IX, 'TIME  STEP  NUMBER= ' ,13/) 
TIME=T*DFLOAT(L) 

WRITEC 6 , 32 )TIME 

32  FORMATC 5X, ’ELAPSED  TIME= F10 . 4 SECONDS '/ ) 
DO  71  1=1, N 

DO  71  J=1,N 

71  THE(I. J)=l . 0-V(I, J) 

IF( N . GT . 1 1 )GO  TO  58 
DO  59  J=1 , N 

WRITEC  6,82) (THECI , J ) , 1=1 , N ) 

82  FORMAT ( 11(2X,F8.6)) 

59  CONTINUE 
GO  TO  54 
58  CONTINUE 

DO  57  J=1,N 

WRITEC 6, 56) (THECI, J),I=1, 11) 

56  FORMATC 1 1 ( 2X , F8 . 6 ) ) 

WRITEC 6, 56 ) (THECI, J), 1=12, N) 

57  CONTINUE 
54  CONTINUE 

ITMAX=0 
20  CONTINUE 
RETURN 
END 


i 
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SOURCE . LAPLAC 

WRITTEN  BY  R.F.  HANDSCHUH 

*****  PROGRAM  #3  ***** 

THIS  PROGRAM  IS  TO  BE  USED  TO  SOLVE  THE  LAPLACE'S  EQUATION 
USING  THE  EXPONENTIAL  FINITE  DIFFERENCE  METHOD. 


IMPLICIT  REAL*8(A-H,0-Z) 

REAL*8  V( 25 , 25 ) 

INPUT  PROGRAM  DATA 

WIIITEC  6 ,15) 

15  FORMAT(  IX, 'NUMBER  OF  NODES=N  I3V) 

READC  9, 10  )N 

10  FORMAT (13) 

WRITEC  6,12) 

12  FORMATC IX, 'NUMBER  OF  TIME  SUB  INTERVALS2  NS  13') 
READC  9 , 1 3 )NS 

13  FORMATC 13) 

WRITEC  6 , 16 ) 

16  FORMATC IX, 'TOTAL  NUMBER  OF  TIME  STEPS2  NTOT  13') 

•  READC  9,21 )NTOT 

21  FORMATC 13) 

WRITEC 6, 24) 

24  FORMATC IX, 'INPUT  TIME  /  LENGTH  SQUARED  F5.3') 

READC  9 , 25 )TSI 

25  FORMATCF5 . 3) 

WRITEC 6, 22) 

22  FORMATC IX, 'TOTAL  TIME  OF  ONE  TIME  STEP2  T  F7.6') 

READC  9 , 23 )T 

23  FORMATC  F7 . 6 ) 

WRITE (6, 26) 

26  FORMATC IX, 'NUMBER  OF  TIME  STEPS  BEFORE  PRINTING2  13’) 
READC  9,27 )IPR 

27  FORMATC 13 ) 

WRITEC 6. 31) 

31  FORMATC IX, 'INPUT  ACCURACY  DESIRED2  F7.6') 

READC  9 , 32 )ET1 

32  FORMATC F7 . 6  ) 

N1 =N- 1 

INITIALIZE  THE  BOUNDRY  CONDITIONS 

DO  30  1=2, N1 
DO  30  J=2,N1 
30  V(I,J)=100. 

C 

I  =  N 

DO  50  J=1,N 
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50  V(I,J)=0.0 
J=N 

DO  51  1=1, N 

51  V(I,J)=0.0 


J=1 

DO  53  1=2, HI 
53  V(I,J)=0.0 


1=1 

DO  52  J=1 , N 
52  VCI,J)=100.0 


CALL  EXP  FIN  DIF  FOR  LAPLACE  EQUATION 

CALL  EFDIHC ( N , NS , NTOT , TSI , V , T , IPR , ET1 ) 

STOP 

END 

SUBROUTINE  EFDIHC 

SUBROUTINE  EFDIHC ( N , NS , NTOT , TSI , V , T , IPR , ET 1 ) 

IMPLICIT  REAL*8(A-H,0~Z) 

REALMS  VT(25,25),V(25,25),l1(25,25)rP(25.25),THE(25,25) 
TS=TSI/DFLOAT ( NS+ 1 ) 

N1=N-1 

NS1=NS+1 

BEGIN  MAIN  STEP  INCREMENT 

DO  20  L=1 , NTOT 
ETA=0 . 0 

ZERO  THE  SUM  OF  THE  DRIVE  NUMBERS 

DO  15  J=1,N 
DO  15  1=1, N 
15  P(I,J)=0. 

SAVE  THE  DEPENDENT  VALUES  FROM  THE  LAST  TOTAL  TIME  STEP 

DO  10  J=1,N 
DO  10  1=1, N 
DELV=V(I,J)-VT(I,J) 

ETA=ETA+DABS( DELV ) 

10  VT(I,J)=V(I,J) 

IFCL.LE.88)  GO  TO  107 
IFCETA.LE.ETDGO  TO  100 
107  CONTINUE 


.*1 


SUB  TIME  INTERVAL 
DO  30  K=1 , NS1 

CALCULATE  THE  DRIVE  NUMBERS 

DO  41  J=2,N1 
JM1=J-1 
JP1=J+1 
DO  40  1=2, N1 
IM1=I-1 

Mci7j)=(VT(IPl,J)+VT(IMl,J)+VT(I,JPl)+VT(I,JMl)-4.*VT(I,J)) 

M(I,J)=M(I,J)/VT(I,J) 

CONTINUE 

CALCULATE  THE  DEPENDENT  VARIABLE  ON  THE  SUB-TIME  INTERVAL 

DO  50  11=2, N1 
DO  50  J1=2,N1 

VT(Il,Jl)  =  VT(Il,jn*DEXP(TS*M(Il,Jl)) 

SUM  THE  DRIVE  NUMBERS 

DO  60  1=1, N1 
DO  60  J=1,N1 
P(I,«7)=P(I»J)+M(I#«J) 

CONTINUE 

CALCULATE  THE  DEPENDENT  VARIABLE  ON  THE  NEXT  TOTAL  STEP 

DO  70  J=1,N 
DO  70  1=1, N 

VCI, J)=V(I, J)*DEXP(TS#P(I, J) ) 

CONTINUE 
ITMAX=ITMAX+ 1 

PRINT  THE  RESULTS 

IF  (UMAX .  LT .  IPR  )GO  TO  20 
WRITE( 6,5) 

FORMATC/) 

WRITE( 6 , 51 )L 

FORMATC IX, 'TIME  STEP  NUMBER=  ,13/) 

TIME=T*DFLOAT(L) 

DO  71  1=1, N 
DO  71  J  =  1 , N 
THEC I, J )=V( I , J) 

IFCN.GT. 11 )GO  TO  58 
DO  59  J=1 , N 

WRITE ( 6 ,82 ) (THEC I , J) ,1=1 ,N) 

FORMATC  1K2X.F8 .4) ) 

CONTINUE 
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GO  TO  54 
58  CONTINUE 

DO  57  J=1,N 

WRITE ( 6 , 56  )  ( THE (I#J)»I=1»11) 

56  FORMAT (11(2X»F8.4)) 

WRITE(6 » 56 ) (THECI, J) , 1=12 ,N) 

57  CONTINUE 

54  CONTINUE  ) 

ITMAX=0  . 

20  CONTINUE  ) 

GO  TO  101  : 

100  CONTINUE  i 

WRITEC6 , 103)L  “ 

103  FORMAT C 2X ' #*###  CONVERGED  RESULT  «*#**', 5X , ’ITERATION* ' ,14// )  I 

DO  102  J=1 ,N 

WRITE  ( 6 ,56)(V(I,J),I=1,N) 

102  CONTINUE  f 

101  CONTINUE  i 

RETURN  ' 

END  5 


a 


a 

5 
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WRITTEN  BY  R.F.  HANDSCHUH 
*****  PROGRAM  #4  ***** 
SOURCE. EFDVAR 


THIS  PROGRAM  IS  FOR  THE  SOLUTION  OF  THE  DIFFUSION  EQUATION 
WITH  VARYING  THERMAL  CONDUCTIVITY.  THE  METHOD  USED  IS  THE 
EXPONENTIAL  FINITE  DIFFERENCE  METHOD. 


IMPLICIT  REAL*8(A-H,0-Z) 

REAL*8  V(100),KR 

INPUT  PROGRAM  DATA 

WRITEC  6 ,15) 

15  FORMAT ( IX, 'NUMBER  OF  NODES=N  13'/) 

READ( 9 , 1 0  )N 

10  FORMAT (13) 

WRITEC  6 ,12) 

12  FORMATC IX, 'NUMBER  OF  TIME  SUB  INTERVALS=  NS  13') 

READ ( 6 , 1 3 )NS 

13  FORMATC 13) 

WRITEC 6, 16) 

16  FORMATC IX, 'TOTAL  NUMBER  OF  TIME  STEPS=  NTOT  13’) 

READC  9,21 )NTOT 

21  FORMATC 13) 

WRITEC 6, 24) 

24  FORMATC IX, 'INPUT  TIME*THERMAL  DIFFUSIVITY  /  LENGTH  SQUARED  F5.3') 
READC  9 , 25 )TSI 

25  FORMATC F5. 3) 

WRITEC 6, 22) 

22  FORMATC IX, 'TOTAL  TIME  OF  ONE  TIME  STEP=  T  F7.5') 

READC  6 , 23  )T 

23  FORMAT (F7 .5) 

WRITEC 6, 14) 

14  FORMATC IX, 'INPUT  REFERENCE  THERMAL  CONDUCTIVITY=F5 . 4  ’  ) 

READC  9,17 )KR 

17  FORMATC F5. 4) 

WRITEC 6, 26) 

26  FORMATC IX, 'INPUT  THERMAL  CONDUCTIVITY  SLOPE  VALUE=  F5 . 4 ’ ) 

READC  9,27 )BETA 

27  FORMATC F5 . 4) 

INITIALIZE  THE  BOUNDRY  CONDITIONS 

VC  1 ) =0  . 

VC N) =0  . 

Ml  =N- 1 

DO  30  1=2, N1 
30  V(I)=100. 
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CALL  EXP  FIN  DIF  FOR  INFINITE  HEAT  TRANSFER  COEFFICENT 

CALL  EFDIHC ( N , NS , NTOT , TSI , V , T , KR , BETA ) 

STOP 

END 

SUBROUTINE  EFDIHC 

SUBROUTINE  EFDIHC ( N , NS , NTOT , TSI , V , T , KR , BETA ) 

FOR  INFINITE  HEAT  TRANSFER  COEFFICENT 


IMPLICIT  REAL*8(A-H,0-Z) 

REAL*8  VT(100),V(100),M(100),P(100) 

REAL*8  KR , K ( 1 0  0 ) , KO ( 1 0  0 ) , THE ( 1 0  0 ) , KAPPA (100) 
TS=TSI/DFLOAT ( NS+ 1 ) 

N1=N-1 

NS1=NS+1 

PRINT  HEADING 

WRITE( 6,5) 

WRITEC6 ,100) 

1 0  0  FORMAT  dx.'*****  SOURCE . EFDVAR  ***** » / ) 
BEGIN  MAIN  TIME  STEP  LOOP 
DO  20  LS1 > NTOT 


ZERO  THE  SUM  OF  THE  DRIVE  NUMBERS 

DO  15  1=1, N 
15  P(I)=0. 


SET  VARIABLES  EQUAL  TO  THE  LAST  STEPS  VALUES 


DO  10  1=1, N 

10  VT(I)=V(I) 

DO  11  1=1, N 

11  K0(I)=KR+BETA*KR*V(I) 

SUB  TIME  INTERVAL 

DO  30  KK=1 , NS1 

DO  35  1=1, N 

K(I ) =KR+BETA*KR*VT( I ) 

35  THE(I)=VT(I)+BETA*VT(I)*VT(I)/2.0 

CALCULATE  THE  DRIVE  NUMBERS 

DO  40  1=2, N1 
IM1=I- 1 
IP1=I+1 
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Md)=VT(IPl)+VT(IMl)-2.*VT(I) 

Md)=M(I)+BETA*<VTdPl)**2.+VTdMl)**2.-2.*VT(I)**2.  )/2. 

40  Md)=M(I)/THE(I) 

DO  50  11=2 , N1 

50  KAPPA(I1 )=THE(I1 )#DEXP(TS*( 1 . +BETA*VT(I1 ) )*MCI1 ) ) 

CALCULATE  THE  DEPENDENT  VARIABLE  ON  THE  SUB-TIME  INTERVAL 
DO  55  11=2, N1 

VT(I1  )  =  (-l .  +SQRTC 1 .  +2  .*KAPPAdl  )*BETA)  )/BETA 

SUM  THE  DRIVE  NUMBERS 

DO  60  1=2, N1 
60  Pd)=Pd)+I1d) 

30  CONTINUE 
WRITE ( 6 , 5 ) 

5  FORMAT ( / ) 

WRITE (6 , 31 )L 

31  FORM AT (IX, 'TIME  STEP  NUMBER= ' , 13/ ) 

TlME=T*DFLOAT(L) 

WRITE( 6 , 32) TIME 

32  FORMAT C5X, 'ELAPSED  TIME=' ,F10 ,4, 'SECONDS'/) 

CALCULATE  THE  NEXT  TOTAL  STEP  DEPENDENT  VARIABLES  AND  PRINT  RESULTS 
DO  70  1=1, N 

THE(I)=Vd)+BETA*Vd)#Vd)/2 . 0 

KAPPA (I ) =THE (I ) *DEXP ( TS* ( 1 . + BETA* V (I ) ) *P Cl ) ) 

70  V(I)=(-1. +SQRTC 1 . +2 . *KAPPA(I)*BETA) )/BETA 

WRITE (6,81)(V(I),I=1,11) 

81  FORMAT ( IX, 11(F8.5,2X) ) 

20  CONTINUE 
RETURN 
END 


non  non  ooo  nnnnn.-innnnnnnn 
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WRITTEN  BY  R.F.  HANDSCHUH 

SOURCE. EXPVAR 

*****  PROGRAM  #5  ***** 

THIS  PROGRAM  IS  FOR  THE  SOLUTION  OF  THE  DIFFUSION  EQUATION 
WITH  VARYING  THERMAL  CONDUCTIVITY.  THE  METHOD  USED  IS  THE 
EXPLICIT  FINITE  DIFFERENCE  METHOD. 


IMPLICIT  REAL*8CA-H,0-Z) 
REAL*8  VCIOOKKR 

INPUT  PROGRAM  DATA 


WRITE ( 6 , 15) 

FORMATC IX,  'NUMBER  OF  NODES=N  I3V) 

READ( 9 , 1 0 )N 
FORMAT(I3) 

WRITE( 6,16) 

FORMAT( IX, 'TOTAL  NUMBER  OF  TIME  STEPS=  NTOT  13’) 

READC 9,21 )NTOT 
FORMATCI3) 

WRITE( 6,24) 

FORMAT( IX, 'INPUT  TIME*THERMAL  DIFFUSIVITY  /  LENGTH  SQUARED  F5.3') 
READC  9 , 25 )TSI 
FORMATCF5.3) 

WRITEC  6,22) 

FORMATC IX, 'TOTAL  TIME  OF  ONE  TIME  STEP=  T  F7 . 5 ’ ) 

READC  6 , 23 )T 
FORMATC F7. 5) 

WRITEC  6,14) 

FORMATC IX, ’INPUT  REFERENCE  THERMAL  CONDUCTIVITY=F5 . 4 ’ ) 

READC  9,17 )KR 
FORMATCF5.4) 

WRITEC  6 , 26 ) 

FORMATC IX, 'INPUT  THERMAL  CONDUCTIVITY  SLOPE  VALUE=  F5 . 4 ' ) 

READC  9,27 )BETA 
FORMAT (F5. 4) 


INITIALIZE  THE  BOUNDRY  CONDITIONS 

VC  1  )  =  0  . 

VCN)=0 . 

N 1  =  H  - 1 

DO  30  1-2 , N1 
V(I)=100 . 


CALL  EXP  FIN  DIF  FOR  INFINITE  HEAT  TRANSFER  COEFFICENT 


onn  ooo  non  non  nno 


CALL  EFDIHC ( N , NTOT , TSI , V , T , KR , BETA ) 

STOP 

END  | 

SUBROUTINE  EFDIHC  [ 

SUBROUTINE  EFDIHC!N, NTOT, TSI, V.T.KR, BETA)  | 

< 

FOR  INFINITE  HEAT  TRANSFER  COEFFICENT  * 

* 

IMPLICIT  REAL*8!A-H,0-Z)  ; 

REAL*8  VT!  100),V(100),M!100),P!100) 

REAL*8  KR  J 

N1=N-1  « 

BEGIN  TIME  STEP  LOOP  f 

DO  20  L=1 ,NTOT  j 

CALCULATE  DEPENDENT  VARIABLE  | 

DO  40  1=2, N1 

IM1 =1-1  f 

IP1=I+1 

VT(I)=V!I)+TSI*! ( 1 . +BETA*V!I) )*!V!IP1 )+V(IMl )-2.*V!I)))  i 

DERSQR=DABS!V(IP1 )-V!IMl ) )  ' 

IF  ( DEP.SQR .  LE .  0 . 0  )  GO  TO  4  0 

VT!I)=VT!I)+TSI*! !BETA*!DERSQR)**2. )/4. )  ! 

40  CONTINUE  < 

i 

PRINT  RESULTS  J 

WRITE( 6 , 5 )  J 

5  FORMAT!/)  < 

WRITE! 6 , 31 )L  i 

31  FORMAT! IX, 'TIME  STEP  NUMBERS ,13/)  \ 

TIME=T*DFLOAT ( L ) 

WRITE! 6 , 32 )TIME 

32  FORMAT !5X, 'ELAPSED  TIME=' ,F10. 4, 'SECONDS'/) 

DO  70  1=1, N 

70  V(I)=VT(I) 

WRITE! 6,81)!V(I),I=1,N) 

81  FORMAT! 1X,11(F9.S,2X))  i 

20  CONTINUE  ! 

RETURN  ' 

END  : 


i 


i 


fo-MU  Ml 


UNCLASSIFIED 


AN  EXPONENTIAL  FINITE  DIFFERENCE  TECHNIOUC  FOR  SOL VI NO  2/ 
PARTIAL  DIFFERENTI.  .  <U>  GAERTNER  <H  U>  RESEARCH  INC 
NORMALK  CT  R  F  HANDSCHUH  JUN  9?  NASA-E-2344 
USAAVSCOH-TR-87-C-1S  F/O  12/2  NL 
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SOURCE .  IMPVAR 


WRITTEN  BY  R.F.  HANDSCHUH 
*****  PROGRAM  #6  ***** 

THIS  ROUTINE  IS  TO  BE  USED  FOR  COMPARISION  TO  EXPONENTIAL 
FINITE  DIFFERENCE  ALGORITHM.  THIS  ROUTINE  WILL  USE  THE 
IMPLICIT  ROUTINE  TO  SOLVE  FOR  THE  TEMPERATURE  FIELD  USING 
THE  TRI-DIAGONAL  MATRIX  ALGORITHM. 


IMPLICIT  REAL*8(A-H,0-Z) 

REALMS  V( 100 ) ,KP 

INPUT  PROGRAM  DATA 

WRITE ( 6,15) 

15  FORMATC IX, ’NUMBER  OF  NOD£S=N  13’/) 

READC 9, 10  IN 

10  FORMAT(I3 ) 

WRITE ( 6 , 16 ) 

16  FORMAT C IX, ’TOTAL  NUMBER  OF  TIME  STEPS=  NTOT  13’) 

READC  9,21 )NTOT 

21  FORMATCI3) 

WRITEC  6,24) 

24  FORMAT ( IX, ’INPUT  TIME#THERMAL  DIFFUSIVITY  /  LENGTH  SQUARED  F5.3’) 
READC  9 , 25 )TSI 

25  FORMAT C  F5 . 3 ) 

WRITEC 6, 22) 

22  FORMAT C IX, ’TOTAL  TIME  OF  ONE  TIME  STEP=  T  F7.5’) 

READC  6 , 23 )TS 

23  FORMAT C  F7 . 5 ) 

WRITEC 6, 14) 

14  FORMATC IX, 'INPUT  REFERENCE  THERMAL  CONDUCTIVITY=F5 . 4 ’  ) 

READC  9 . 1 7 )KR 

17  FORMATC F5. 4) 

WRITEC 6, 26) 

26  FORMATC IX, ’INPUT  THERMAL  CONDUCTIVITY  SLOPE  VALUE=  F5.4’) 

READC  9,27 )BETA 

27  FORMATC F5. 4) 

INITIALIZE  BOUNDRY  CONDITIONS 

VC1)=0. 

VCN)  =  0  . 

N1=N-1 

DO  30  1=2, N1 
30  VCI)=100. 

CALL  IMPLICIT  ROUTINE 


nnno'f!  non  n  non  o  non  non  non  non 
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CALL  IMPL  <  N , NTOT , TSI , V , TS , KR , BETA ) 

STOP 

END 

SUBROUTINE  IMPL 

SUBROUTINE  IMPL ( N , NTOT , OMEGA , T , TS , KR , BETA ) 

FOR  INFINITE  HEAT  TRANSFER  COEFFICENT  AND  VARYING  THERMAL  CONDUCTIVITY 
IMPLICIT  REAL*8 ( A-H  » 0-Z ) 

REAL*#  A(101),B(101),C(101),D(101),KAPC101),GAM(101),TC101) 

REAL*8  TO(lOl) 

REAL*8  KR 

PRINT  HEADING 

WRITE (6  *  331 ) 

331  FORMATC IX, '**#*«  SOURCE . IMPVAR  *#*##*/) 

RHO=l . 0 
CP=1 . 0 
N1=N-1 

BEGIN  TIME  STEP  LOOP 
DO  20  L=1 . NTOT 

CALCULATE  THOMAS  ALGORITHM  VARIABLES  AND  THOSE  THAT  ARE  A  FUNCTION 
OF  TEMPERATURE 

DO  21  1=1, N 
21  D(I)=T(I) 

DO  200  1=1, N 
200  TO(I)=T(I) 

DO  25  1=2, N1 
KAP(I)=1.0+BETA*TCI) 

25  GAM(I)=BETA*(T(I+l)-T(I-l))/(4. ) 

DO  30  1=2, N1 

A(I)=-OMEGA*(KAP(I)-GAM(I) ) 

B(I)  =  (l.+2. *OMEGA*KAP ( I ) ) 

30  C(I)=-OMEGA*(KAP(I)+GAM(I)) 

CALL  TRI-DIAGONAL-MATRIX  ALGORITHM 

CALL  TRIDAG(2,N1,A,B,C,D,T) 

0  CONTINUE 

PRINT  THE  RESULTS 

WRITEC 6,5) 

5  FORMATC /) 


nnn  nnn  nnnnnnnn 
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WRITE(6,31)L 

31  FORMAT (IX,' TIME  STEP  NUMBERS ,13/) 

TIME=TS*DFLOAT ( L ) 

WRITE C  6 > 32 ) TIME 

32  F0RMAT(5X, ’ELAPSED  TIME= ’, FI  0 . 4 , 'SECONDS’ / ) 

WRITE ( 6 ,81)(T(I) , 1=1 , N) 

81  FORMAT ( 1X,11(F9.5,2X) ) 

20  CONTINUE 
RETURN 
END 

SUBROUTINE  TRIDAG 

THIS  ROUTINE  IS  FOR  THE  SOLUTION  OF  THE  THOMAS  ALGORITHM 

THIS  ROUTINE  WAS  TAKEN  FROM  THE  BOOK  APPLIED  NUMERICAL  METHODS 
BY  CARNAHAN ,  LUTHER,  AND  WILKES. 

SUBROUTINE  TRIDAG (IF,L,A,B,C,D,V) 

IMPLICIT  REAL*8 ( A~H , O-Z ) 

REAL*8  A(101),B(101),C(101),D(101),V(101) ,BETA( 101 ) , GAMMA ( 101 ) 

COMPUTE  INTERMEDIATE  ARRAYS  BETA  AND  GAMMA 

BETA(IF)=BCIF) 

GAMMA ( IF )=D( IF) /BETA (IF) 

IFP1=IF+1 
DO  1  I=IFP1 , L 

BETA(I)=B(I)-A(I)*C<I-1)/BETACI-1) 

1  GAMMA(I)=(D(I)-A(I)*GAMMA(I-1 ) )/BETA(I) 

COMPUTE  FINAL  SOLUTION  VECTOR  V 

V(L)=GAMMA(L) 

LAST=L-IF 
DO  2  K= 1 . LAST 
I  =  L-K 

2  V( I) =GAMMA( I)-C(I )*V(I+1 )/BETA( I) 

RETURN 

END 


oonn  ono  onnonnonoonnnn 
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SOURCE . COUE 

WRITTEN  BY  R.F.  HANDSCHUH 
*****  PROGRAM  *7  ***** 

THIS  PROGRAM  IS  TO  BE  USED  TO  DEMONSTRATE  THE  USE  OF  THE 
EXPONENTIAL  FINITE  DIFFERENCE  METHOD  ON  THE  DEVELOPING  TEMPERATURE 
FIELD  IN  A  LAMINAR  COUETTE  FLOW 


IMPLICIT  REAL«8 (A-H.O-Z) 

REAL*8  V(100) 

INPUT  PROGRAM  DATA 

WRITE ( 6,15) 

15  FORMAT( IX, 'NUMBER  OF  NODES=N  13' /) 

READC 9 , 1 0  )N 

10  FORMAT( 13 ) 

WRITE( 6,12) 

12  FORMAT( IX, 'NUMBER  OF  SUB  INTERVALS3  NS  13’) 

READ( 6,13 )NS 

13  FORMAT  C 1 3 ) 

WRITEC  6,161 

16  FORMAT ( IX, 'TOTAL  NUMBER  OF  POSITION  STEPS3  NTOT  13') 

READC  9,21 )NTOT 

21  FORMATCI3) 

WRITEC 6, 24) 

24  FORMATC IX, 'INPUT  (POSITION  STEP#KIN.  VISC)/(LENGTH  SQUARED)  F5.3’) 
READC  9 , 25 )TSI 

25  FORMAT CF5. 3) 

WRITEC 6. 22) 

22  FORMATC IX, 'TOTAL  DISTANCE  OF  ONE  STEP3  T  F5.3’) 

READC  6 , 23 )T 

23  FORMATCF5 .3) 

WRITEC 6. 14) 

14  FORMATC IX, 'INPUT  (MAX  WIDTH)/(FLOW  VEL)  F5.4') 

READC  6 , 1 7 )SP 

17  FORMATC  F5 . 4 ) 

WRITEC 6 . 31 ) 

31  FORMATC IX. 'INPUT  INTERVAL  FOR  PRINTING  RESULTS3  13’) 

READC  6 , 32 )IPR 

32  FORMATC 13) 


INITIALIZE  THE  BOUNDRY  CONDITIONS 

VC1)  =  0  . 

DO  30  1  =  2, N 
30  V(I)=100. 
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CALL  EXPONENTIAL  FINITE  DIFFERENCE  FOR  COUETTE  FLOW 

CALL  EFFL ( N , NS , NTOT , TSI , V , T , SP , IPR ) 

STOP 

END 

SUBROUTINE  EFFL 

SUBROUTINE  EFFL ( N , NS , NTOT , TSI , V , T , SP , IPR ) 

FOR  COUETTE  LAMINAR  FLOW 
IMPLICIT  REAL#8 (A-H.O-Z) 

REAL#8  VT(100) ,V(100) ,H(100),P(100),FL(20) 
TS=TSI/DFL0AT(NS+1 ) 

N1=N-1 


PRINT  HEADING 
WRITE C  6 ,  92 ) 

92  FORMAT! IX, 'SOLUTION  FOR  DEVELOPING  TEMPERATURE  FIELD  IN' 

* ' LAMINAR  COUETTE  FLOW'//) 

DY=l./DFLOAT(N-l) 

CALCULATE  PARAMETER  THAT  VARIES  WITH  POSITION  IN  THE  FLOW 
DO  115  1=2, N 

L 1 5  FL ( I ) =SP/ ( DY*DFLOAT ( I- 1 ) ) 

NS1=NS+1 

BEGIN  TOTAL  POSITION  STEP  LOOP 

DO  20  L=1 , NTOT 
IT=IT+ 1 

ZERO  THE  SUM  OF  DRIVE  NUMBERS 

DO  15  1=1, N 
15  PCI)=0. 

SET  TEMPORARY  VALUES  EQUAL  TO  THE  LAST  POSITION  STEP  VALUE 

DO  10  1=1, N 
10  VT(I)=V(I) 

SUB  POSITION  INTERVAL 

DO  30  K= 1 , NS1 

CALCULATE  THE  DRIVE  NUMBERS 


DO  AO  1=2, N1 
IM1 =1- 1 


oor  nno  non 


IP1=I+1 

40  M(I>  =  (2.*VTCI)-VT(IM)-VT(IP1))/VT(I) 

CALCULATE  TEMPORARY  DEPENDENT  VARIABLE  ON  SUB-INTERVAL 
DO  50  11=2, N1 

50  VT(I1)=VT(I1) #DEXP ( -TS#M(I1 )#FL(I1 ) ) 

DO  60  1=2, N1 
60  P(I)=P(I)*M(I) 

30  CONTINUE 

CALCULATE  THE  COMPLETE  STEP  DEPENDENT  VARIABLES 
DO  70  1=1, N 

V(I)=V(I)#DEXP(-TS#P(I)#FL(I) ) 

70  CONTINUE 

IF ( IT . LT . IPR ) GO  TO  20 
IT=0 

PRINT  THE  RESULTS 

WRITE ( 6,5) 

5  FORMATC/) 

WRITE ( 6 , 31 )L 

31  FORMATC IX, 'POSITON  STEP  NUMBER =', 13) 

DIST=T#DFLOAT ( L ) 

WRITE (6,32) DIST 

32  FORMATC5X, 'LOCATION  DIST= F10 . 4 . 'METERS ' ) 

WRITE ( 6 , 82)(V(I)»I=1,N) 

82  FORMATC 1X,11(F9.5,2X)) 

20  CONTINUE 
RETURN 
END 


nnn  nnn  nnonnnnnonnoonn 


WRITTEN  BY  R.  F.  HANDSCHUH 
SOURCE. EX3D 

*****  PROGRAM  #8  ***** 


THIS  ROUTINE  IS  FOR  FINDING  THE  TEMPERATURE  AT  A  GIVEN  LOCATION 
AND  TIME  FOR  A  3-DIMENSIONAL  SOLID. 

THIS  PROGRAM  CALCULATES  THE  EXACT  TEMP  AS  FOUND  IN  THE  BOOK 
"TRANSPORT  PHENOMENA"  BY  BIRD,  STEWART,  AND  LIGHTFOOT. 


IMPLICIT  REAL*8! A-H.O-Z) 

REAL*8  PI,PI2,PI3 
PI=3. 1415926 

INPUT  PROGRAM  DATA 

URITE( 1,31) 

31  FORMAT( IX, ’NUMBER  OF  NODES  PER  COORDINATE  DIRECTION=I3 ’ ) 
READ( 5 , 32  INODES 

32  FORMATCI3) 

WRITE! 1,1) 

1  FORMAT! IX. 'INPUT  INITIAL  TEMPERATURE  FOR  THE  SOLID  F6  3’) 
READ! 5 , 2 )T0 

2  FORMAT !F6. 3) 

WRITE! 1,3) 

3  FORMAT! IX, ’SURFACE  TEMPERATURE  FOR  TIME  >  0  SECONDS  F6.3’) 
READ! 5 , 4 )T1 

4  FORMAT !F6. 3) 

WRITE! 1.5) 

5  FORMAT! IX, ’INPUT  THERMAL  DIFFUSIVITY  F5.3’) 

READ! 5, 6) ALFA 

6  FORMAT !F5. 3) 

ASSUMING  3  DIM.  CUBE 

DATA  A, B, C/1. 0,1. 0,1.0/ 

WRITE! 1.9) 

9  FORMAT! IX. ’INPUT  THE  NUMBER  OF  TIMES  THROUGH  SUMMATION  13’) 
READ! 5 , 1 1 )N1 

11  FORMAT! 13 ) 

WRITE! 1,12) 

12  FORMAT! IX. 'INPUT  TIME  TO  BE  EVALUATED  AT  F5.3’) 

READ! 5,14 )TIME 

14  FORMAT! F5. 3) 

WRITE! 6 , 20 )T0 ,T1 

20  FORMAT! IX. ’TEMP  INITIAL= ’ , El  5 . 8 , 2X , ’TEMP  SURF= ’ , El 5 . 8/ ) 
WRITE! 6 , 3  0 ) ALFA , TIME 

30  FORMAT! IX, ’THER  DIFF= ’ , El 5 . 8 , 2X , • TIME= ’ , F7 . 5 ) 


non  o  nno  non 
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START  SUMMATION  ROUTINE 

PI3=PI**3. 

PI2*PI **2. 

CALCULATE  THE  TEMPERATURE  ALONG  THE  DIAGONAL 

DEL=A/(DFLOAT( NODES- 1) ) 

DO  100  NODE2 1, NODES 
X=DEL*DFLOAT(NODE-l ) 

Y=X 

sum=o.o 

DO  10  K=1 ,N1 
RP=DFLOAT(K-l) 

RP1=RP+ . 5 

COSP =DCOS ( RP 1 #PI*Z/C ) 

DO  10  N=1 , N1 
RN=DFLOAT(N-l ) 

RN1=RN+ . 5 

COSN2DCOS(RNl#PI#Y/B) 

DO  10  M=1.N1 
RM=DFLOAT(M-l) 

RM1=RM+ . 5 

COSM=DCOS ( RM1 *PI#X/A ) 

J=M+N+K-3 
Jl=J/2 
J2=J1*2 
VAL=-1 . 0 

IF( J . EQ . J2 ) VAL=1 . 0 

GAM2(RM1#*2)/(A*A)+(RN1**2)/(B*B)+(RP1**2)/(C*C> 
EXPLIM= ( -GAM*PI2*ALFA*TIME ) 

IFCEXPLIM.LT. -100. )GO  TO  999 
EP=DEXP(EXPLIM) 

GO  TO  998 
999  CONTINUE 
EP=0 . 0 

998  CONTINUE 

SUM=SUM+(VAL/(RM1*RN1*RP1*PI3) )*EP#COSM*COSN*COSP 
10  CONTINUE 

TNEW=Tl+8 . #(T0~T1 )*SUM 

PRINT  THE  RESULTS 

WRITE C 6 , 16 )X. Y>Z 

16  FORMATC IX , ' X= ' , F5 . 3 , 2X» ’ Y= ' . F5 . 3 , 2X , *Z= ’ , F5 . 3 ) 
IFCTMEW . LT . IE- 1 0 )TNEW=0 . 0 
WRITE ( 6 , 1 5 )TNEW 

1 5  FORMAT ( 1 X . ’ TEMPERATURE2 ' , El  5 . 8 ) 

100  CONTINUE 
STOP 
END 


v.vlv.  ^V-V.V.V. 
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SOURCE. EFD 3D 

WRITTEN  BY  R.F.  HANDSCHUH 
*****  PROGRAM  #9  ***** 

THIS  PROGRAM  IS  FOR  3-DIMENSIONAL  CARTESIAN  COORDINATES 

UNSTEADY  STATE  HEAT  TRANSFER  IN  A  CUBE.  THE  METHOD  OF  SOLUTION  IS  THE 

EXPONENTIAL  FINITE  DIFFERENCE  ALGORYTHM.  THIS  PARTICULAR  PROGRAM  IS 

FOR  INFINITE  HEAT  TRANSFER  COEFFICENT  AT  THE  EXPOSED  SURFACES 

AT  X2Y2Z=1.Q  FOR  THE  THREE  SURFACES  WHERE  X,Y,Z  EQUAL  0.0 

ARE  TO  BE  CONSIDERED  AS  PERFECTLY  INSULATED. 


IMPLICIT  REAL*8C A-H.O-Z) 

REAL*8  VC  25 , 25 , 25 ) 

INPUT  PROGRAM  DATA 

WRITE ( 6 .15) 

15  FORMATC IX, 'NUMBER  OF  NODES=N  13’/) 

READC 9 , 1 0 )N 

10  FORMAT(I3) 

WRITE( 6,12) 

12  FORMATC IX, 'NUMBER  OF  TIME  SUB  INTERVALS2  NS  13') 

READC  9 , 1 3 )NS 

13  FORMATC 13) 

WRITE C  6,16) 

16  FORMATC IX, 'TOTAL  NUMBER  OF  TIME  STEPS2  NTOT  13’) 

READC  9,21 )NTOT 

1  FORMATC 13) 

WRITEC  6,24) 

24  FORMATC IX, ’INPUT  TIME*THERMAL  DIFFUSIVITY  /  LENGTH  SQUARED  F5.3’) 
READC  9 , 25 )TSI 

25  FORMATC F5. 3) 

WRITEC 6, 22) 

22  FORMATC IX, 'TOTAL  TIME  OF  ONE  TIME  STEP2  T  F6.4') 

READC  9 , 23 )T 

23  FORMAT CF6 . 4 ) 

WRITEC 6, 26) 

26  FORMATC IX, ’NUMBER  OF  TIME  STEPS  BEFORE  PRINTING2  13’) 

READC  9 , 27 )IPR 

27  FORMATC 13) 

WRITEC 6 , 250 )N, NS, NTOT 

250  FORMATC IX,’#  OF  NODES2 ’, 13 , 2X t  OF  SUB-TIME-INT2 ' , 13 , 2X, 

*'#  OF  TIME  STEPS2 ’,13) 

WRITEC  6,251 )TSI,T 

251  FORMATC IX, ' CTIME*THER  DIFF) /LENGTH  SQUARED2 ’ F5 . 3 , 2X , 

* ' TIME  STEP  LENGTH2' ,F6 .4/) 

N12N-1 

INITIALIZE  BOUNDRY  CONDITIONS 


oo  non  nno  non  noo  non  nnn 
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DO  30  1=1. HI 
DO  30  J=1,N1 
DO  30  K=1,N1 
30  V(I,J,K)=1.0 

C 

I=N 

DO  50  J=1 . H 
DO  50  K=1,H 

50  V(I,J,K)=0.0 
J=N 

DO  51  1=1, N 
DO  51  K=1.N 

51  V(I,J,K)=0.0 

DO  52  1=1, N 
DO  52  J=1 , N 

52  V(I,J,K)=0.0 

CALL  EXP  FIN  DIF  FOR  INFINITE  HEAT  TRANSFER  COEFFICENT 

CALL  EFDIHC C  N , NS , NTOT , TSI , V , T , IPR ) 

STOP 

END 

SUBROUTINE  EFDIHC 

SUBROUTINE  EFDIHC(N , NS , NTOT , TSI , V , T , IPR ) 

FOR  INFINITE  HEAT  TRANSFER  COEFFICENT 
IMPLICIT  REAL#8(A-H,0-Z) 

REAL*8  VT(25,25,25),V(25,25,25),M(25,25,25),P(25,25,25) 
REAL*8  Ml. M2 
TS=TSI/DFLOAT ( NS+ 1 ) 

N1=N-1 

NS1=NS*1 

PRINT  HEADING 

WRITE( 6 ,200) 

200  FORMATC1X, ***********  RESULTS  FROM  EFD3D  #**#####»* » // ) 
START  TOTAL  TIME  STEP  LOOP 
DO  20  L=1 , NTOT 

ZERO  THE  SUM  OF  THE  SUB-INTERVAL  DRIVE  NUMBERS 

DO  15  K=1,N 
DO  15  J= 1 , N 
DO  15  1=1, N 
15  P( I , J , K ) =0 . 

SET  SUB-INTERVAL  VALUES  EQUAL  TO  THE  LAST  TIME  STEP  VALUES 


A 


DO  10  K=1 ,  N 
DO  10  J=1,N 
DO  10  1=1, N 
VT(I,J,K)=V(I,J,K) 

SUB  TIME  INTERVAL 

CALCULATE  THE  DRIVE  NUMBERS  WHICH  IS  DEPENDENT  ON  LOCATION  IN  THE  CUBE 

DO  30  KS=1 , NS1 

DO  42  K=2,N1 

KM1  =K- 1 

KP1=K+1 

DO  41  J=2,N1 

JM1=J-1 

JP1=J+1 

DO  40  1=2, N1 

IM1  =  I-1 

IP1=I+1 

M1=VT(IP1, J,K)*VT(IM1 , J,K)+VT(I,JP1 ,K)+VT(I, JMl.K) 

M2=M1+VT(I, J.KPl )+VT(I, J,KMl)-6 .*VT(I,J,K) 

IF(VT(I,J,K) . LE .0.0)MCI,J,K)=0. 

IF( VT( I , J, K) . LE . 0 . 0 )GO  TO  40 
MCI. J,K)=M2/VT(I,J,K) 

CONTINUE 

CONTINUE 

CONTINUE 

INSULATED  BOUNDRY  ALONG  X-AXIS 

J  =  1 
K=1 

KP1=K+1 
JP1=J+1 
DO  48  1=2, N1 
IP1=I+1 
IM1=I-1 

M1=VT(IP1 , J,K)+VT(IM1 ,J,K)+2.#VTCI, JP1 ,K) +2 . *VT(I, J, KP1 ) 
IF(VT(I,J,K).LE.0.0)M(I,J,K)*0. 

IF(VT(I, J,K) .LE. 0 . 0)GO  TO  48 

MCI, J,K)=(Ml-6 . #VT(I. J,K) )/VT(I,J,K) 

CONTINUE 

INSULATED  BOUNDRY  ALONG  Y-AXIS 
1=1 

IP1=I+1 
K=  1 

KP1 =K+ 1 
DO  43  J=2,N1 
JP1=J+1 
JM1 :  J-  1 

Ml =2 . *VT ( IP1 , J , K ) +VTC I , JP1 , K ) + VT( I , JM1 , K ) +2 . * VT( I , J , KP3 ) 


nnn  ono  ooo  cion 


A 


IFCVTCI,J,K).LE.0.0)rKI,J,K)=0. 

IF(VT(I,J,K) .LE.O .0 )G0  TO  43 
M(I,J,K)=(Ml-6 .#VT(I,J,K) )/VT(I,J,K) 

43  CONTINUE 

INSULATED  BOUNDRY  ALONG  Z-AXIS 
J=1 

JP1=J+1 

1=1 

IP1=I+1 
DO  44  K=2,N1 
KM1=K-1 
KP1 =K+ 1 

ril=2.*VT(IPl,J,K)+2.*VT(I, JP1 ,K)+Vf(I, J,KP1)+VT(I, J.KM1) 

IF( VT( I >  J » K ) . LE.0.0)M(I,J,K)=0. 

IF( VT( I , J » K ) .LE.O . 0 )GO  TO  44 
M(I,J,K)  =  (t11-6. #VT(I, J , K) )/VT(I, J,K) 

44  CONTINUE 

INSULATED  FACE  AT  Z=0 
K=  1 

KP1=K+1 

DO  45  1=2, N1 

IP1=I+1 

IM1 =1“ 1 

DO  45  J=2.N1 

JP1=J*1 

JM1=J-1 

m=VT(IPl,J,K>+VTCIMl,J,K)+VT(I,UPl,K)+VT(I,Jt11,K)  +  2.*VTCI,J,KPl) 
IF(VT(I,J,K) .LE. 0.0)M(I,J,K)=0. 

IF( VT( I . J, K ) . LE . 0 . 0 )GO  TO  45 
tt(I,J,K)=ail-6.*VT(I,J,K))/VTCI,J,K> 

45  CONTINUE 

AT  THE  FACE  WHERE  Y=0 
J=1 

JP1=J+1 

DO  46  1=2, N1 

IP1=I+1 

I«1=I-1 

DO  46  K=2,N1 

KP1=K+1 

Kttl =K- 1 

t11=VT(IPl, J,K)+VT(IM1, J,K)+2.*VT(I, JP1 ,K ) +VTCI , J , KP1 ) +VT(  I ,  J ,  KM1 ) 
IF(VT(I,  J,K)  .LE.  0. 0)WI,  J,K)  =  0. 

IF(VT(I,J,K> .LE.O. 0)GO  TO  46 

M(I, J,K)=(Ml-6 .*VT(I, J,K))/VT(I,J,K> 

46  CONTINUE 

AT  THE  FACE  WHERE  X=0 
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! 

» 

i 

s 

i 

j 


i=i 

ipi=i+i 

DO  47  J=2,N1 
JP1=J+1 

jni=j-i 

DO  47  K=2 » N1 

KP1=K+1 

KM1=K-1 

Ml=2.*VT(IPl,J,K)+VTCI,JPl,K)+VT(r,  JM1 ,K)+VTCI, J,KP1 )+VT(I,J,KMl ) 
IFC  VTC I ,  J  ,K )  .LE.0.0)MCI,J,K)=0. 

IF(VTCI,J,K).LE.O.O)GO  TO  47 
mi,j,K)  =  (m-6.*VT(i,j,ta)/VT(i,j,K) 

47  CONTINUE 
C 

C  CORNER  AT  ORIGIN 
C 

Ml  =  2.*VTCl,2,n*2.*VTC2,l,l)+2.*VTCl,l,2)-6.*VT(  1,1,1) 
flf  l.l*l)*ni/VT(l,l,l) 

C 

C 

C 

C  CALCULATE  THE  SUB-INTERVAL  DEPENDENT  VARIABLES 

C 

DO  50  11=1, N 
DO  50  Jl=l ,N 
DO  50  K1=1.N 

IF(MCI1,J1,K1) .LT.-50. ) VTCI1 , J1 ,K1 )=0 . 0 
.  IFCMCI1,J1,K1).LT.-5Q)G0  TO  50 

VTCI1,J1,K1)=VTCI1,J1,K1)*DEXP(TS*M(I1,J1,K1)) 

50  CONTINUE 
C 

C  SUM  THE  DRIVE  NUMBERS 

C 

DO  60  1=1, N 
DO  60  J=1,N 
DO  60  K= 1 , N 

60  PCI, J,K)=PCI,J,K)+M(I,J,K) 

30  CONTINUE 
C 

C  CALCULATE  THE  NEXT  COMPLETE  TIME  STEP  DEPENDENT  VARIABLES 
C 

DO  70  K=1,N 
DO  70  J=1,N 
DO  70  1=1, N 

IFCPCI. J,K) .LT.-50. Wl,  J,K)=0. 0 
IFCPCI.J.K) .LT.-50. JGO  TO  70 
V(I,J,K)=VCI, J,K)*DEXP(TS*P(I,J,tO) 

70  CONTINUE 

ITMAX=ITMAX+ 1 

IK( ITMAX . LT . IPR )GO  TO  20 

WRITEC  6 , 5 ) 

5  FORMATC / ) 

WRITEC  6 , 31 )L 

31  FORMATC IX, 'TIME  STEP  NUMBER =’ ,13/) 


onn 
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TIME=T#DFLOAT(L) 

WRITEC 6  >  32 )TIME 

32  F0RHATC5X, 'ELAPSED  TIflEs ’ » FI 0 . 4 , ' SECONDS ' / ) 
PRINT  OUT  THE  DIAGONAL  RESULTS 


WRITEC 6, 82 )( V( I ,1,1), 1=1, N) 
82  FORMAT (11C2X.F8.6)) 

ITMAXr 0 
20  CONTINUE 
RETURN 
END 
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SOURCE. EXPL3D 

WRITTEN  BY  R.F.  HANDSCHUH 

*****  PROGRAM  #10  ***** 


THIS  PROGRAM  IS  FOR  3-DIMENSIONAL  CARTESIAN  COORDINATES 

UNSTEADY  STATE  HEAT  TRANSFER  IN  A  CUBE.  THE  METHOD  OF  SOLUTION  IS  THE 

PURE  EXPLICIT  METHOD.  THIS  PARTICULAR  PROGRAM  IS 

FOR  INFINITE  HEAT  TRANSFER  COEFFICENT  AT  THE  EXPOSED  SURFACES 

AT  X=Y=Z=1 . 0  FOR  THE  THREE  SURFACES  WHERE  X.Y.Z  EQUAL  0.0 

IS  TO  BE  CONSIDERED  AS  PERFECTLY  INSULATED. 


IMPLICIT  REAL*8( A-H.O-Z) 

REAL*8  V(25,25,25) 

INPUT  PROGRAM  DATA 

NUMBER  OF  NODES  =  N 

READ( 5 ,  1 0  )N 
10  FORMAT( 13 ) 

TOTAL  NUMBER  OF  TIME  STEPS  =  NTOT 

READC 5,21 )NTOT 
21  FORMAT (14) 

TIME#THERMAL  DIFFUSIVITY/LENGTH  SQUARED  =  TSI 

READ( 5 . 25 )TSI 
25  FORMAT(F5.3) 

TOTAL  TIME  OF  ONE  TIME  STEP  =  T 

READC  5 , 23 )T 
23  FORMATCF6 .4) 

NUMBER  OF  STEPS  BEFORE  PRINTING  THE  RESULTS  =  IPR 

READ( 5 , 27 )IPR 
27  FORMATC 13 ) 

WRITEC6 ,252) 

252  FORMATC IX, • a##*###*####*###*#  PURE  EXLICIT  FINITE  DIFFERENCE  '  - 
^'METHOD  *****************'///) 

WRITEC  6 , 250 )N , NS , NTOT 

250  FORMATC  IX,*#  OF  NODES3  '  .  13 , 2X ,  '  #  OF  SUB-TIME-INT=  M3 , 2X , 

#'#  OF  TIME  STEPS3 ' > 14 ) 

WRITEC 6, 251 )TSI,T 

251  FORMATC IX, ' (TIMEMTHER  DIFF)/LENGTH  SQUARED2 ' F5 . 3 » 2X , 

* ' TIME  STEP  LENGTH3  * ,F6 . 4/  ) 

N1=N-1 


onn  nnoono  0000000  000  o  000 
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INITIALIZE  BOUNDS!  CONDITIONS 

DO  30  1=1, N1 
DO  30  J=1 ,N1 
DO  30  K=1 ,N1 
30  V(I. J,K)=1 . 0 

I=N 

DO  50  J=1 , N 
DO  50  K=1,N 

50  V(I.J,K)=0.0 
J-N 

DO  51  1=1, N 
DO  51  K=1,N 

51  V(I, J,K)=0 . 0 
K-N 

DO  52  1=1, N 
DO  52  J=1,N 

52  V(I, J,K)=0 . 0 

CALL  PURE  EXPLICIT  FINITE  DIFFERENCE  FOR  INFINITE  HEAT  TRANSFER  COEFFICENT 

CALL  PURE( N , NTOT , TSI , V , T , IPS ) 

STOP 

END 

SUBROUTINE  PURE 

SUBROUTINE  PURE ( N , NTOT , TSI , V , T , IPR ) 

PURE  EXPLICIT  FINITE  DIFFERENCE  METHOD 
FOR  INFINITE  HEAT  TRANSFER  COEFFICENT 

IMPLICIT  REAL*8(A-H.O-Z) 

REALMS  V(25.25,25),VT<25,25,25) 

REAL*8  Ml, M2 
N1=N-1 

START  TIME  STEP  LOOP 
DO  20  L= 1 , NTOT 

SAVE  VALUES  FROM  THE  LAST  TIME  STEP 

DO  39  1=1, N 
DO  39  J=1,N 
DO  39  K=1,N 

39  VT(I,J,K)=V(I,J,K) 

CALCULATE  THE  FIELD  VARIABLE  USING  THE  EXPLICIT  FINITE  DIFFERENCE  METHOD 

DO  42  K=2,N1 
KM1=K-1 


non  non  non 
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KP1=K+1 

DO  41  J=2 ,N1 

JM1=J-1 

JP1=J+1 

DO  40  1=2, HI 

IM1=I-1 

IP1=I+1 

M1=VT(IP1,J,K)+VT(IP11,J,K)+VT(I,JP1,K)+VT(I,JP11,K) 
M2=M1+VT(I , J,KP1)+VT(I#J, KM1 ) -6 . *VT( I, J , K ) 
V(I,J,K)=V(I,J,K) +TSI#M2 

40  CONTINUE 

41  CONTINUE 

42  CONTINUE 

INSULATED  BOUNDRY  ALONG  X-AXIS 

J  =  1 

K  =  1 

KP1=K+1 
JP1=J+1 
DO  48  1=2, N1 
IP1=I+1 
IM1=I-1 

M1=VT(IP1,J,K)+VT(IM1,J,K)+2.*VT(I,JP1,K)+2.#VT(I,J,KP1) 
M2=Ml-6 .*VT(I,J,K) 

V(I.J,K)=V<I,J,K)+TSI*ri2 
48  CONTINUE 

INSULATED  BOUNDRY  ALONG  Y-AXIS 
1=1 

IP1=I+1 

K=1 

KP1=K+1 
DO  43  J=2,N1 
JP1=J+1 
JM1=J-1 

Ml=2.*VT(IPl,J.K)+VT(I,JPl,K)+VT(I,jm,K)+2.*VTCI,J,KPl> 
M2=Ml-6 . *VT(I, J,K) 

V(I,J,K)=V(I,J,K) +TSI*M2 

43  CONTINUE 

INSULATED  BOUNDRY  ALONG  Z-AXIS 


J=1 

JP1=«J  +  1 
1  =  1 

IP1=I+1 
DO  44  K=2,N1 
KM1 =K- 1 
KP1 =K+ 1 

Ml =2 .*VT(IP1,J,K)+2.«VT(I,JP1,K)+VT(I,J,KP1)+VT(I,J,KM1) 
M2=Ml-6 .*VT(I, J.K) 

V(I,J,K)=V(I,J,K) +TSI*t12 


V.Y.V.  V'w.V.V 
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44  CONTINUE 
INSULATED  FACE  AT  Z=0 
K=1 

KP1=K+1 

DO  45  1=2, N1 

IP1=I+1 

IM1=I-1 

DO  45  J=2,N1 

JP1=J+1 

JM1=J-1 

M1=VT(IP1 , J,K)+VT(IM1 , J,K) +VT( I, JP1 ,  K)+VT(I ,  JPll , K)+2 . #VT(I, J,KP1 ) 
M2=M1-6.*VT(I,J,K) 

V(I.J,K)=V(I,J,K)+TSI#M2 

45  CONTINUE 

AT  THE  FACE  UHERE  Y=0 
J=1 

JP1=J+1 

DO  46  1=2, N1 

IP1=I+1 

IM1=I-1 

DO  46  K=2,N1 

KP1=K+1 

KM1=K-1 

M1=VT(IP1 'J,K)+VT(IMl,J,K)+2.*VT(InJPl,K)+VT(IiJrKPl)+VT(I,J,Kni ) 
P12=Ml-6  .  #VT(I, J ,K) 

V(I,J,K)=V(I,J,K)+TSI*M2 

46  CONTINUE 

AT  THE  FACE  UHERE  X=0 
1=1 

IP1=I+1 

DO  47  J=2,N1 

JP1=J+1 

JM1=J-1 

DO  47  K=2 ,N1 

KP1=K+1 

KM1=K-1 

m=2.#VT(IPl,J,K)+VT(I,JPl,K)+VT(I, Jt11  ,K)+VT(I.  J,KP1  J+VTCI,  J,KM1 ) 
t12=Ml-6.*VT(I,J,K) 

V(I»J,K)=V(I,J,K)  +TSI*f12 

47  CONTINUE 

CORNER  AT  ORIGIN 

M1=2.»VT(1,2,1)+2.*VT(2,1,1)+2.*VT(1,1,2)-6.*VT( 1,1,1) 
V(l,l,l)=V(l,l,l)+m*TSI 


ITP1AX=ITMAX+1 


« 
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IF ( ITP1AX . LT . IPR ) GO  TO  20 
URITEC  6,5) 

5  FORMATC/) 

URITEC  6 , 31  )L 

31  fORMATClX. 'TIME  STEP  NUMBER®' ,13/) 
TIME=T*DFLOAT(L) 

URITEC 6 , 32 )TIME 

32  FORMATC5X, 'ELAPSED  TIME® ’ ,F10 . 4 , 'SECONDS '/ ) 

PRINT  OUT  THE  DIAGONAL  RESULTS 

URITEC 6, 82 )C VC  I, I, I), 1=1, N) 

82  FORMATC 11C2X.F8.6) ) 

ITMAX- 0 
20  CONTINUE 
RETURN 
END 


non  ooo  nonno.innnnnnnnon 
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SOURCE. DOUGLA 


WRITTEN  BY  R.F.  HANDSCHUH 
*****  PROGRAM  *11  ***** 

THIS  PROGRAM  IS  FOR  3-DIMENSIONAL  CARTESIAN  COORDINATES 

UNSTEADY  STATE  HEAT  TRANSFER  IN  A  CUBE.  THE  METHOD  OF  SOLUTION  IS  THE 

METHOD  OF  DOUGLAS.  THIS  PARTICULAR  PROGRAM  IS 

FOR  INFINITE  HEAT  TRANSFER  COEFFICENT  AT  THE  EXPOSED  SURFACES 

AT  X2Y=Z2 1 . 0  FOR  THE  THREE  SURFACES  WHERE  X.Y.Z  EQUAL  0.0 

IS  TO  BE  CONSIDERED  AS  PERFECTLY  INSULATED. 


IMPLICIT  REAL*8( A-H.O-Z) 

REAL*8  T(25, 25 , 25 ) 

INfUT  PROGRAM  DATA 

WRITEC6.15) 

15  FORMAT C IX, ’NUMBER  OF  NODES 2 N  13' /) 

READC  5  ,  1 0  )N 
10  FORMAT (13) 

WRITE( 6,16) 

16.  FORMAT ( IX. 'TOTAL  NUMBER  OF  TIME  STEPS=  NTOT  13’) 

READ(5,21) NTOT 

21  FORMAT(I3) 

WRITEC 6.24) 

24  FORMATC IX, 'INPUT  TIME*THERMAL  DIFFUSIVITY  /  LENGTH  SQUARED  F5.3’) 
READC  5,25 )TSI 

25  FORMAT( F5 . 3 ) 

WRITEC  6 , 22 ) 

22  FORMAT ( IX, 'TOTAL  TIME  OF  ONE  TIME  STEP=  T  F6.4') 

READC5,23)DT 

23  FORMATC F6. 4) 

WRITEC 6. 26) 

26  /ORMAT ( IX, 'NUMBER  OF  TIME  STEPS  BEFORE  PRINTING2  13’) 

READC  5 . 27 )IPR 

27  FORMATC 13 ) 

WRITEC 6 ,250 )N, NS, NTOT 

250  FORMATC IX, '#  OF  NODES2' ,13, 2X, '*  OF  SUB-TIME-INT2' ,13, 2X, 

*’#  OF  TIME  STEPS2' ,13) 

WRITEC6 ,251 )TSI , DT 

251  FORMATC  IX,  '  (TT.ME*THER  DIFF)/LENGTH  SQUARED2  'F5. 3, 2X, 

* ' TIME  STEP  LENGTH2' ,F6 .4/) 

N1=N-1 


INITIALIZE  BOUNDRY  CONDITIONS 

DO  30  1=1. N1 
DO  30  J=1,N1 
DO  30  K=1,N1 


noon  ooo  non  non  ooo  non 
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30  T(I,J,K)=1.0 

C 

I=M 

DO  50  J=1 » N 
DO  50  K=1,N 

50  T(I, J»K)  =  0 . 0 
J=N 

DO  51  1=1, H 
DO  51  K=1 , H 

51  T(I,J,K)=0.0 

I/-U 

DO  52  1=1, N 
DO  52  J=I  ,N 

52  T(I,J,K>=0.0 
C 

CALL  DG(N,NS.NTOT,TSI,T,DT,IPR) 

:;top 

END 

SUBROUTINE  DG 

SUBROUTINE  DG ( N , NS , NTOT , TSI , T , DT , IPR ) 

FOR  INFINITE  HEAT  TRANSFER  COEFFICENT 
IMPLICIT  REAL*8( A-H.O-Z) 

REALM  U(25,25,25),VC25»25,25),M(25,25,25),T(25,25,25> 

REALMS  TEMP (25) ,A(25),B(25),C(25),D(25) 

REALM  Ml, M2 
N1=N-1 

PRINT  OUT  HEADING 
WRITE( 6,200) 

200  FORMATC IX,  •##########  RESULTS  FROM  METHOD  OF  DOUGLAS  *#####*###'  - 
E// ) 

CALCULATE  COEFFICENTS  FOR  THOMAS  ALGORITHM  (  TRI-DIAGONAL  MATRIX  SOLVER) 

A(  1  )  =  0  .  0 
B( 1 ) = 1 . +TSI 
C( 1 ) = -TSI 
DO  60  1=2, N 
A( I ) =- . 5*TSI 
B( I ) = 1 . +TSI 
60  C(I)=A(I) 

BEGIN  TIME  STEP  LOOP 

DO  20  L=1 , NTOT 

CALCULATE  TEMPORARY  VARIABLES  "U"  -  X  DIRECTION,  ""V"  -  Y  DIRECTION, 

THEN  ACTUAL  FIELD  VARIABLE  "T"  -  Z  DIRECTION. 


nan  nnn 
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DO  120  L00P=1 , 3 

DO  42  K=2 ,  N1 

DO  41  J=2,N1 

DO  40  1=2, N1 

IF ( LOOP . EQ . 2 ) GO  TO  140 

IF ( LOOP . EQ . 3 ) GO  TO  141 

DELX=T(I+1,J,K)+T(I-1,J,K)-2.*TCI,J,K) 

DELY=T(I,d*l,K)+T(I, J- 1 ,K ) -2 . *T( I, J ,K ) 

DELZ=T(I,0,KU)+T(I,J,K-1)-2.*T(I,J,K) 

M( I , J » K ) =T( I , J » K)  +TSI* ( .  5#DELX+DELY +DELZ ) 

GO  TO  40 

140  CONTINUE 

DELU=U(in,J,K)+U(I-l,J,K)-2.«U(I,J,K) 
DELY=T(I,J+1,K)+T(I,J-1,K)-2.#T(I, J,K) 
WI,J,K)=M(I,J,K)*.5#TSI*<DELU-DELY) 

GO  TO  40 

141  CONTINUE 

DELV=V(I,J+1,K)+V(I,J-1,K)-2.*V(I,J,K) 
DELZ=T(I,d,K+l)+T(I,J,K-l)-2.*T(I, J,K) 

H(I,J,K)=H(I,  J,K)+.5*TSI#(DELV-DELZ) 

40  CONTINUE 

41  CONTINUE 

42  CONTINUE 

INSULATED  BOUNDRY  ALONG  X-AXIS 

J  =  1 
K=1 

KP1=K+1 
JP1=J+1 
DO  48  1=2, N1 
IP1-I+1 
1111  =  1-1 

IF ( LOOP .  EQ .  2 )GO  TO  148 
IF ( LOOP . EQ . 3 ) GO  TO  248 

M1=.5#(T(I+1,J,K)+T(I-1,J,K)-2.*T(I,J,K» 

M2=2.#T(I, JP1,K)+2.#T(I.J,KP1)-4.#T(I,J,K) 
mi,J,K)=T(I,J,K)+TSI*(ni+M2) 

GO  TO  48 
148  CONTINUE 

DELU=U(I+1 ,J,K)+U(I-1,J»K)-2.*U(I»J,K) 
DELY=2.*CT(I.J+1,K)-T(I,J,K)) 

M(I,J,K)=M(I. J,K)+.5#TSI#(DELU-DELY) 

GO  TO  48 
248  CONTINUE 

DELV=2 . *V(I,J+1,K)-2.*V(I,J,K) 

M(I»J»K)=M(I.J»K)+. 5#TSI#( DELV-2 .*(T(I,J,KP1)-T(I,J,K))) 
48  CONTINUE 

INSULATED  BOUNDRY  ALONG  Y-AXIS 


1  =  1 

IP1*I+1 
K=  1 
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KP1=K+1 
DO  43  J=2,N1 
JP1=J+1 
jru*J-i 

IF(LOOP.EQ. 2) GO  TO  143 
IF ( LOOP . EQ . 3 ) GO  TO  243 

m=T(IPl,J,K)+T(I,JPl,K)+T(I,Jm,K)+2.*TCI,J,KPl)-5.*TCI,J,IO 

M(I.J.K)=T(I»J,K)+TSI*m 

GO  TO  43 

143  CONTINUE 

DELU=2.*(U(IP1,J,K)-U(I,J,K» 

DELY=T(I.JP1,K)+T(I, JM1 ,K)-2 . *T(I, J,K) 

M( I , J , K ) =t1( I . 0 , K) ♦ . 5*TSI* ( DELU-DELY ) 

GO  TO  43 

243  CONTINUE 

DELV=V  (I,\J+1,K)+V(I,J-1,K)-2.*V(I,J,K) 

M(I, J,K)=M(I,J,K)+.5*TSI»(DELV-2.*(T(I,J,KP1)-T(I,0M<))) 

43  CONTINUE 

INSULATED  BOUNDRY  ALONG  Z-AXIS 
J=1 

JP1=J-H 

1=1 

IP1=I+1 
DO  44  K=2,N1 
KM1=K-1 
KP1=K+1 

IF (LOOP . EQ . 2 )GO  TO  144 
IFCLOOP.EQ. 3) GO  TO  244 

m=T(IPl,J,K)  +  2.*T(I,JPl,K)+T<I,J,KPl)+T(I,J,Kt11)-5.*T(I,J,K) 

MI,J»K)>T(X,J,K)+TSX*m 

GO  TO  44 

144  CONTINUE 
m=nci.J,K) 

M2=M1+.5*TSI*(2.*(U(IP1,J,K)-U(I,<J,K))-2.#(T(I,JP1,K)-T(I,J,K))) 

n(I,J,K)=M2 

GO  TO  44 

244  CONTINUE 

DELZ=T(I,J.K+1)+T(I,J»K-1)-2.*T(I,J,K) 

M(I,v7.K)=t1(I.J.K)+. 5*TSI*( 2.#(V(I»JP1,K)-V(I,J,K)) -DELZ ) 

44  CONTINUE 

5  INSULATED  FACE  AT  Z=0 

r* 

K  =  1 

KP1=K+1 
DO  45  1=2  »  N1 
IP1=I+-1 

im=i-i 

DO  45  J=2,N1 

UP1=J+1 

JM1=J-1 

IFCLOOP.EQ. 2)GO  TO  145 
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IF ( LOOP . EQ . 3 ) GO  TO  245 

m=.5*(T(IPl, J,K)+T(IM1,J,K)-2.#T(I,J,K)) 
M2=T(I,JP1,K)+T(I,JM,K)  +  2.*T(I,J,KP1)-4.*TCI,J,K) 
n(I,J,K)=T(I,J,K)+TSI*(m+!12) 

GO  TO  45 

145  CONTINUE 

DELU=U(I+1,J,K)+U(I-1,J,K)-2.#U(I,J,K) 
DELY=T(I,J+1,K)+T(I» J-l ,K)-2 .*T(I, J.K) 
M(I,J,K)=ri(I,J,K)+.5*TSI#(DELU-DELY) 

GO  TO  45 

245  CONTINUE 

DELV=V(I,J+1,K)+V(I,J-1,K)-2.#V(I,J,K) 

MI,J,K)=M(I,J,K)+.5*TSI#(DELV-2.*(T(I,J,KPn-T(I,J,K))) 

45  CONTINUE 

AT  THE  FACE  WHERE  Y=0 
J=1 

JP1=J+1 

DO  46  1=2, N1 

IP1=I+1 

im=i-i 

DO  46  K=2,N1 

KP1 =K+1 

KM1=K~1 

IF( LOOP . EQ . 2 )GO  TO  146 
IFC LOOP . EQ . 3 )GO  TO  246 

m=.5*(T(IPl,J,K)+T(im,J,K)-2.#T(I,J,K)) 
M2=2.*T(I,JP1,K)+T(I,J,KP1 l+TCI, J,KM1 )-4.#T(I,J,K) 
M(I,J,K)=T(I,J.K)+TSI*(Pll*M2) 

GO  TO  46 

146  CONTINUE 

DELU=U(I+1,J.K)+U(I-1,J,K)-2.*U(I,J,K) 

M(I,J,K)=H(I»J,K)+. 5*TSI*(DELU~2 .*(T(I,JP1»K)-T(I,J»K))) 
GO  TO  46 

246  CONTINUE 

DELZ=T(I,J,K+1)+T(I,J,K-1)-2.*T(I,J,K) 

M(I,J,K)=M(I,J,K)+. 5#TSI#( 2 .*(V(I,JP1,K)~V(I,J,K) )-DELZ) 

46  CONTINUE 

AT  THE  FACE  WHERE  X=0 


f)f)n 


GO  TO  4 7 
147  CONTINUE 

DELY=T(I, J+1,K)+Td, J-l ,K) -2 . #T( I , J ,K ) 

n(I,J.K)=M(I. J,K)+.5*TSI*(2.#(U(IPl,J,K)-U(I,J,IO)-DELY) 

GO  TO  47 
247  CONTINUE 

DELV=V(I,J+1 ,K)+Vd,J-l,K)-2.#Vd,J,K) 

DELZ=Td.J,K*l)*Td. J,K-l)-2.*Td, J,K) 

Md,  J,K)=Md,  J,K)+.5«TSI#(DELV-DELZ) 

47  CONTINUE 

CORNER  AT  ORIGIN 

IF C  LOOP . EQ . 2 )  GO  TO  151 
IF ( LOOP . EQ . 3 )GO  TO  251 

M1=T(2,1,1)+2.«T(1»2.1)*2.*T(1,1,2)-5.*T(1,1,1> 
M(1,1,1)=T(1.1.1 )+TSI*Ml 
GO  TO  51 
151  CONTINUE 

M1=MU.1,1> 

M2=Ml  +  .5*TSI*(2.#<U(2,l,l)-UU,l,l)>-2.*(TU,2,l)-Td,l,l))) 
mi,i,n=n2 
GO  TO  51 
251  CONTINUE 

mi=mu.i  ,i) 

M2=M1 ♦ . 5*TSI#( 2 .#( V( 1 ,2, l)-V(l,l,l))-2.*(Td»l,2)_Td,l»l))) 
M(  1 , 1 , 1 )  =M2 
51  CONTINUE 
C 

c 

IF( LOOP . EQ . 2 )GO  TO  130 
IFC  LOOP . EQ . 3 )GO  TO  230 
C 

DO  70  K= 1 , N1 
DO  70  J=1,N1 
DO  30  1=1, N1 
30  D( I) =M( I. J,K) 

C 

C  CALL  TRI  DIAGONAL  MATRIX  ALGORITHM 
C 

CALL  TP.IDAGU  ,N1  , A , B ,C , D . TEMP ) 

C 

DO  430  1=1, HI 
430  UCI, J,K)=TEMP(I) 

C 

70  CONTINUE 
GO  TO  330 
130  CONTINUE 

DO  71  K= 1 , N 1 
DO  71  I -  1 , N1 
DO  72  J=1 , N 1 
72  D(J)=Md,J,K) 
c 

CALL  TRIDAG ( 1 ,N1 , A , B , C , D, TEMP) 


DO  472  J=1,N1 
472  V!I,J,K)=TEMP!J) 

71  CONTINUE 
GO  TO  330 
230  CONTINUE 

DO  73  1=1, HI 
DO  73  J=1,N1 
DO  74  K  =  1 , N 1 
>4  D!K)=M(I. J,K) 

CALL  TRIDAG( 1 ,N1 , A , B , C , D, TEMP) 

DO  474  K= 1 . N1 
474  T(I, J,K)=TEMP!K) 

73  CONTINUE 
330  CONTINUE 
120  CONTINUE 

ITMAX=ITMAX+ 1 

IF(  UMAX  .  LT  .  IPR  )GO  TO  20 

WRITE( 6,5) 

5  FORMAT!/) 

WRITE! 6 , 91 )L 

91  FORMAT! IX, 'TIME  STEP  NUMBER* 1 ,13/) 

TIME=DT*DFLOAT!L) 

WRITE! 6 , 32 ) TIME 

32  FORMAT! 5X, ’ELAPSED  TIME= *, FI 0 . 4 ,’ SECONDS V ) 

PRINT  OUT  THE  DIAGONAL  RESULTS 

WRITE! 6,82)(T(I,I,I),I=1,N) 

82  FORMAT! 11!2X,F8.6) ) 

ITMAX=0 
20  CONTINUE 
RETURN 
END 

SUBROUTINE  TRIDAG 

SUBROUTINE  TRIDAG! IF , L , A , B , C , D, V ) 

IMPLICIT  REAL*8! A-H.O-Z) 

REAL*8  A!25),B!25) ,C!25),D(25) , V! 25 ) , BETA! 25 ) , GAMMA! 25 ) 

COMPUTE  INTERMEDIATE  ARRAYS  BETTA  AND  GAMMA 

BETA!IF)=B!IF) 

GAMMA! IF )=D(IF)/BETA!IF) 

IFP1=IF+1 
DO  1  I=IFP1,L 

BETA !I)=B!l)-A!I)*C!I-l ) /BETA! I- 1 ) 

1  GAMMA! I )=!D!I)-A!l) #GAMMA! 1-1 ) ) /BETA! I ) 
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I 

COMPUTE  FINAL  SOLUTION  VECTOR  V 

V(L)=GAMMA(L) 

LAST=L-IF 
DO  2  K=1 .LAST 
I=L-K 

2  V(I)=GAMMA(I)-C(I)#V(I+1)/BETA(I) 

RETURN 
END 
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WRITTEN  BY  R.F.  HANDSCHUH 

SOURCE . BURGER 

*****  PROGRAM  *12  ***** 

THIS  PROGRAM  IS  FOR  THE  SOLUTION  OF  BURGER'S  EQUATION 
BY  THE  EXPONENTIAL  FINITE  DIFFERENCE  METHOD. 


IMPLICIT  REAL*8CA-H,0-Z) 

REAL*#  VC  100) 

READ  IN  DATA  TO  BE  USED  IN  THE  SOLUTION 
WRITEC  6  >15) 

15  FORMATC IX, 'NUMBER  OF  NODES=N  I3V) 

READC  9 , 1 0 )N 

10  FORMATC 13) 

WRITEC6. 12) 

12  FORMATC IX, 'NUMBER  OF  TIME  SUB  INTERVALS-  NS  13') 

READC 9, 13)NS 

13  FORMATCI3 ) 

WRITEC 6, 16) 

16  FORMATC IX, 'TOTAL  NUMBER  OF  TIME  STEPS3  NTOT  13') 

READC 9,21 )NTOT 

21  FORMATC 13) 

WP.ITEC  6,29) 

24  FORMATC IX, 'INPUT  TIME  /  LENGTH  SQUARED  F5.3') 

READC  9 , 25 )TSI 

25  FORMATCF5. 3) 

WRITEC 6, 26) 

2 6  FORMAT  C 1 X . ' INPUT  KINEMATIC  VISCOSITY3  F5 . 3 ’ ) 

READC  9,27 )RNU 

27  FORMAT CF5. 3) 

WRITEC 6. 22) 

22  FORMATC IX, 'TOTAL  TIME  OF  ONE  TIME  STEP3  T  F5.3') 

READC  9 , 23 )T 

23  FORMATC F5. 3) 

WRITEC 6, 14) 

14  FORMATC IX, 'INPUT  NUMBUR  OF  STEPS  BETWEEN  PRINTS=I3') 
READC  9,17 )IPR 

17  FORMATC 13) 

DATA  FOR  INITIAL  AND  BOUNDRY  CONDITIONS 

VC  1 )=0  . 

VC  N  )  =  1  . 

N1=N-1 

DO  30  1=2 , N1 
30  VCI)=1. 

CALL  EXPONETIAL  FINITE  DIFFERENCE  FOR  BURGER'S  EQUATION 


c 

CALL  BURG ( N , NS , NTOT , TSI . V , T , RNU , IPR ) 

STOP 

END 

C 

C  SUBROUTINE  BURG 
C 

SUBROUTINE  BURG ( N , NS , NTOT , TSI , V , T , RNU , IPR ) 

C 

c 

IMPLICIT  REAL*8 ( A-H , O-Z ) 

REAL*8  VT(100),V(100),M(100)#P(100) ,THE( 100 ) 
l'S=TSI/DFLOAT(NS+l ) 

DX  =  1 . 0/DFLOAT(N~1) 

N1 =N~  1 
NS  1=  HSU 
WRITE ( 6,900) 

900  FORMAT!/ '  *****  SOLUTION  FOR  BURGER  EQUATION 
C 

C  TOTAL  TIME  STEP  LOOP 
C 

DO  20  L=l,NTOT 
ITMAX=ITMAX+1 
C 

C  ZERO  THE  SUM  OF  DRIVE  NUMBERS 
C 

DO  15  1=1, N 
15  P(I)=0. 

C 

C  SET  THE  TEMPOARY  FIELD  VARIABLE  EQUAL  TO  THE  LAST  TIME  STEP  VALUE 
C 

DO  10  1=1, N 
10  VT(I)=V(I) 

C 

C  SUB  TIME  INTERVAL 

C 

DO  30  K=1 , NS1 
C 

C  CALCULATE  THE  SUB-INTERVAL  DRIVE  NUMBERS 

C 

DO  40  1=2, N1 
IM1 =1" 1 
IP1=I+1 

IF( VT( I) . LE . 0 . 0 )GO  TO  40 

M(I)='0.5*DX*(1.-VT<I))*(VTCIP1)-VTCIM1))/VT(I) 
M(I)=M(I)+RNU*(VT(IP1 )+VT(IMl)-2 .*VT(I) )/VT(I) 

40  CONTINUE 

C 

C  CALCULATE  THE  SUB-INTERVAL  DEPENDENT  VARIABLES 

C 

DO  50  11=2, N1 
CHECK=TS*M( II ) 

IK ( CHECK. LE. -50 . )VT(I1 )=0 .0 
IF ( CHECK . LE . -50 . )GO  TO  50 
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VT  < 1 1 ) = VT ( II ) *DEXP ( TS*M( II ) ) 
50  CONTINUE 

SUM  THE  DRIVE  NUttBERS 

DO  60  1=2, N1 
60  P(I)=P(I)+M(I) 

30  CONTINUE 


CALCULATE  THE  DEPENDENT  VARIABLE  AT  THE  NEXT  COMPLETE  STEP 

DO  70  1=1, N 
CHECK=TS#P(I) 

IFCCHECK . LE . -50 . )V(I)=0 . 0 
IF (CHECK . LE . -50 . )GO  TO  70 
V(I)=V(I)*DEXP(TS*P(I) ) 

70  CONTINUE 

OUTPUT  THE  RESULTS 

IF ( ITMAX . LT . IPR ) GO  TO  20 

ITMAX=0 

WRITE ( 6,5) 

5  FORHAT(/) 

WRITE ( 6 , 31 )L 

31  FORMATUX, 'TIME  STEP  NUMBER=’,I3) 

TIME=T*DFLOAT(L) 

WRITE (6, 32) TIME 

32  FORMAT(5X, ’ELAPSED  TIME= ', FI 0 . 4, 'SECONDS' > 

ISTEP=(N-1)/10 

DO  110  1=1, N 
110  THE(I)=1 . 0-V(I) 

DO  80  J=1,ISTEP 

IS=(J-1)*11+1 

IFIN=J*10-H 

WRITE( 6 . 8 1 ) ( THE (I ) , I=IS , IFIN ) 

81  FORMAT ( 1X,11(F8.6,2X)) 

80  CONTINUE 
20  CONTINUE 
RETURN 
END 
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WRITTEN  BY  R.F.  HANDSCHUH 

SOURCE . EXBURG 

*****  PROGRAM  *13  ##*## 

THIS  PROGRAM  IS  FOR  THE  SOLUTION  OF  BURGER'S  EQUATION  USING  AN  EXPLICIT 
TECHNIQUE.  THE  RESULTS  WILL  BE  USED  TO  COMPARE  TO  THE  EXPONENTIAL 
FINITE  DIFFERENCE  TECHNIQUE. 

IMPLICIT  REAL*8 ( A-H » O-Z ) 

REAL*8  V(100) 

INPUT  PROGRAM  DATA 

WRITE(6  ,15) 

FORMAT(  IX, 'NUMBER  OF  NODES=N  I3V) 

READ( 9 . 1 0  )N 
FORMAT( 13 ) 

WRITE ( 6,16) 

FORMAT ( IX, 'TOTAL  NUMBER  OF  TIME  STEPS=  NTOT  13') 

READC  9,21 )NTOT 
FORMAT(I3) 

WRITE ( 6,24) 

FORMAT ( IX, 'INPUT  TIME  /  LENGTH  SQUARED  F5.3') 

READC  9 , 25)TSI 
FORMATCF5 . 3 ) 

WRITEC  6 , 26 ) 

FORMAT ( IX, 'INPUT  KINEMATIC  VISCOSITY=  F5.3') 

READC  9 , 27 )RNU 
FORMAT CF5. 3) 

WRITEC  6,22) 

FORMATC IX, 'TOTAL  TIME  OF  ONE  TIME  STEP=  T  F5.3’) 

READC  9,23)T 
FORMAT CF5. 3) 

WRITEC 6, 14) 

FORMATC IX, 'INPUT  NUMBER  OF  STEPS  BETWEEN  PRINTS=I3') 

READC  9,17 )IPR 
FORMATCI3 ) 

INITIALIZE  THE  BOUNDRY  CONDITIONS 

VC1)  =  1 . 

VC  N)  =  0  . 

N1=N-1 

DO  30  1=2, N1 
VCI)  =  0  . 

CALL  EXPLICIT  FINITE  DIFFERENCE  SOLUTION  FOR  BURGER’S  EQUATION 

CALL  BURG  C  N , NTOT , TSI , V , T , RNU , IPR ) 

STOP 
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END 

SUBROUTINE  BURG 

SUBROUTINE  BURG ( N , NTOT , TSI , V , T , RNU , IPR ) 


IMPLICIT  REAL*B(A-H,0-Z) 

REAL*#  VT(100),V(100),THE<100) 

PRINT  HEADING 

URITEC6 , 999 ) 

999  FORMAT ( IX, *****  EXPLICIT  BURGER  S  EQT  SOLUTION  ****’/) 
DX=1 . O/DFLOATCN-l ) 

N1=N-1 

TIME  STEP  LOOP 

DO  20  L=1 , NTOT 
ITMAX=ITMAXH 

SAVE  LAST  TIME  STEP  VALUES 

DO  10  1=1, N 
10  VT(I)=V(I) 

EVALUATE  EXPLICIT  FINITE  DIFFERENCE  EQUATION 

DO  40  1=2, N1 

IM1=I-1 

IP1=I+1 

V( I) =VT(I)-VT(I)*T*( VTCIPl )-VT(IMl ) )/(2 . *DX) 
V(I)=V(I)+RNU#T*<VT(IP1)+VT(IM1)-2.*VT(I>)/(DX*DX) 

40  CONTINUE 

WRITE  OUT  THE  RESULTS 

IF (ITMAX.LT. IPR) GO  TO  20 

ITMAX=0 

WRITE ( 6 , 5 ) 

5  FORMAT(/) 

WRITEC  6 , 31  )L 

31  FORMAT(lX, 'TIME  STEP  NUMBER*', 13) 

TIME=T*DFLOAT ( L ) 

WRITE( 6 , 32 )TIME 

32  FORMAT(5X, 'ELAPSED  TIME* ’, F10 . 4 , ’SECONDS’ ) 
ISTEP=(N-1)/10 

DO  110  1=1. N 
110  THE(I)=V(I) 

DO  #0  J= 1 , ISTEP 
IS=(J-1)*11+1 
IFIN=J#1 0+1 

WRITE ( 6 , 8 1 ) ( THE (I) ,I=IS, IFIN ) 
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SOURCE. NONBOU 

WRITTEN  BY  R.F.  HANDSCHUH 

#####  PROGRAM  #14  #«»«« 

THIS  PROGRAM  IS  USED  FOR  THE  SOLUTION  OF  THE  BOUNDRY  LAYER  FLOW 
OVER  A  FLAT  PLATE.  THE  DIRECTION  OF  FLOW  IS  IN  THE  X-DIRECTION 
WHICH  IS  USED  AS  THE  MARCHING  DIRECTION  FOR  THE  EXPONENTIAL  FINITE 
DIFFERENCE  ALGORITHM.  THE  THERMAL  AND  VELOCITY  BOUNDRY  LAYERS 
CAN  BE  EXTRACTED  FROM  THE  TEMPERATURE  AND  VELOCITY  FIELDS  FOUND. 


IMPLICIT  REAL*8! A-H.O-Z) 

REAL*8  U!101),T!101),V!101) 

INPUT  PROGRAM  DATA 

WRITE!  6,15) 

15  FORMAT(  IX, 'NUMBER  OF  NODES  IN  Y  DIRC=N  13'/) 

READ( 9 , 1 0  )N 

10  FORMAT! 13) 

WRITE! 6 , 12 ) 

12  FORMAT! IX, 'NUMBER  OF  SUB  INTERVALS3  NS  13') 

READ! 9 , 1 3 )NS 

13  FORMAT!I3) 

WRITE! 6, 16) 

16  FORMAT! IX. 'TOTAL  NUMBER  OF  X-DIR  STEPS*  NTOT  13') 
RE AD ! 9 , 2 1 ) NTOT 

21  FORMAT! 13) 

WRITE! 6 ,24) 

24  FORMAT! IX. 'INPUT  STEP  LENGTH  F5.3’) 

READ! 9 , 25 )DX 

25  FORMAT! F5. 3) 

WRITE! 6, 26) 

26  FORMAT! IX. 'NUMBER  OF  STEPS  BEFORE  PRINTING*  13') 
READ! 9 , 27 )IPR 

27  FORMAT! 13) 

WRITE! 6, 110) 

110  FORMAT! IX, 'INPUT  KINEMATIC  VISCOSITY*  F6 . 4 ' ) 

READ! 9,111)  RNU 

111  FORMAT! F6. 4) 

WRITE! 6,101) 

101  FORMAT! IX. 'INPUT  THERMAL  DIFFUSIVITY  F6 . 4 ' ) 

READ! 9 , 1 02 )RAL 

102  FORMAT! F6. 4) 

WRITE! 6,103) 

103  FORMAT! IX, 'INPUT  YMAX  F5.1’) 

READ! 9, 104)  YMAX 

104  FORMAT! F5 . 1 ) 

WRITE ! 6 , 2  5  0 ) N , NS , NTOT 

250  FORMAT! IX,'#  OF  NODES* ', 13 , 2X, '#  OF  SUB-INT* ’ , 13 , 2X, 
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*'#  OF  TIME  STEPS= ' ,13) 

RN=DFL0AT(N-1 ) 

DY=YMAX/RN 
WRITE ( 6 , 251 )DX, DY 

251  FORMAT ( IX, *  DX= ' , F8 . 4, 2X, '  DY= ' , F8 . 4/ ) 

WRITE ( 6 , 252 )RNU , RAL 

252  FORMAT( IX, ’KINEMATIC  VISCOSITY= ' , F6 . 4 , 'CM*CM/S * , 2X, 
* ' THERMAL  DIFFUSIVITY= ’ , F6 . 4 , ’ CM*CM/S ’ ) 

INITIALIZE  BOUNDRY  CONDITIONS 

DO  30  J=2,N 
U( J)  =  l .  0 
V(J)=0.0 
30  T( J)  =  l .  0 

U(1)=0.0 
TC 1 )=0 . 0 


CALL  NON1 ( N , NS , NTOT , RNU , U , V , T , IPR , DX , DY , RAL ) 

STOP 

END 


SUBROUTINE  NON 1 ( N , NS , NTOT , RNU , U , V , T , IPR , DX , DY , RAL ) 


IMPLICIT  REAL#8 ( A-H , O-Z ) 

REAL*8  U(lOl) , MU (101) ,V(101) , T ( 1 0 1 ) , MT (101) 

REAL*8  PU(101),PT(101),UT(101) ,TT( 1 0 1 ) , VT( 101 ) 

REAL*8  THE (101,1000, 3), UT 1(101) 

PRINT  HEADING 

WRITE ( 6 , 5 ) 

WRITEC  6 , 222 ) 

222  FORMATC1X, 'a###*####*#  SOURCE . NONBOU  a**#**###**#'//) 

WRITE( 6,223) 

223  FORMAT(/. IX, ’SOLUTION  FOR  BOUNDRY  LAYER  FLOW  PAST  A  FLAT  PLATE’/) 
DY2=DY*DY 

TS=DX#RAL/(DY2*DFLOAT(NS+l ) ) 

TS1 =DX*RNU/ ( DY2*DFLOAT( NS+ 1 ) ) 

DEL=DX/DFLOATCNS+l ) 

HI =N- 1 

NS1=NS+1 

NSTEP=CN-1)/10 

BEGIN  TOTAL  INTERVAL  LOOP  FOR  L=1  TO  NTOT  STEPS 
DO  20  L= 1 , NTOT 

ZERO  THE  SUM  OF  DRIVE  NUMBERS  FOR  THE  NEXT  SET  OF  SUB-POSITION  INTERVALS 


A 


c 

DO  15  1=1, N 
PU(I)=0.0 
15  PT(I)=0 . 0 

C 

C  SAVE  LAST  POSITION  STEP  VALUES  TOR  TEMPORARY  VARIABLE  CALCULATIONS 
C  ON  SUB-POSITION  INTERVAL 

C 

DO  10  1=1, N 
VT(I)=V(I) 

UT(I)=U(I) 

10  TT(I)=T(I) 

C 

C  SUB  -  POSITION  INTERVAL 

C 

DO  30  K=1 , NS1 
C 

C  CALCULATE  TEMPERATURE  FIELD  DRIVE  NUMBER 

C 

DO  41  J=2,N1 

JM1=J-1 

JP1=J+1 

MT(J)=-VT(J)*(TT(JP1)-TT<JM1))*DY/(2.*RAL*UT(J)*TTCJ)) 

MT( J ) =MT( J)+ (TT( JP1 )+TT( JM1 )-2 . *TT( J) )/(UT( J)#TT( J) ) 

41  CONTINUE 

C 

C  CALCULATE  X  -  DIRECTION  VELOCITY  DRIVE  NUMBER 
C 

DO  141  J=2,N1 

JM1=J-1 

JP1=J+1 

MU(J)=-.5*VT(J)*DY*(UT(JPl)-UT(JMl))/<RNU*UTCJ)#UTC>n) 

MUC J)=MU( J)+ ( UT( JP1 ) +UT( JM1 ) -2 . *UT( J) )/( UT( Jl^UTC J) ) 

141  CONTINUE 

C 

C  CALCULATE  TEMPERATURE.  X-DIRECTION  VELOCITY,  AND  Y-DIRECTION  VELOCITY 
C  ON  THE  SUB -POSITION  INTERVAL 
C 

DO  50  11=2, N1 

50  TT(I1 )=TT(I1 )*DEXP(TS*MT(I1 ) ) 

C 

DO  51  1=2, N1 
UT1(I)=UT(I) 

51  UT(I)=UT(I)*DEXP(TS1#MU(I) ) 

C 

DO  65  J=2.N1 
JM1=J-1 

VTCJ)=VT(JM1)-.5*(DY/DEL)*CUT(J)-UT1(J)+UT(JM1)-UT1(JM1)) 

65  CONTINUE 

C 

C  SUM  THE  DRIVE  NUMBERS 

C 

DO  60  J=2,N1 
PU ( J ) =PU ( J ) +MU ( J ) 
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jp  wwi-  *<  A*  *jiff«RMWMi/  v'juuuwi 
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60  PTCJ)=PTCJ)+MTCJ) 

30  CONTINUE 

CALCULATE  THE  NEXT  TOTAL  POSITION  STEP  VALUES  OF  VELOCITIES  AND  TEMPERATURE 

DO  70  J=1 ,  N 
UT(J)=U(J) 

UC J)=UC J)*DEXPCTS1#PUC J) ) 

TCJ)=TCJ)#DEXPCTS#PTCJ)) 

70  CONTINUE 

DO  75  J=2 , N1 
JM1=J-1 

7  5  VCJ)=VCJM1)~. 5*CDY/DX) #CUC J ) -UTC J )  +  UC  JM1 ) -UTC JM1 ) ) 

ITMAX=ITMAX+1 

SAVE  THE  VALUES  FOUND  IN  3-DIMENSIONAL  ARRAY  "THE" 

DO  76  J=1  ,N 
THECJ.L, 1)=UCJ) 

THECJ.L, 2)=VCJ) 

76  THECJ.L, 3)=TCJ) 

IF ( ITMAX . LT . IPR ) GO  TO  20 

WRITE  OUT  THE  RESULTS  AT  THE  REQUESTED  INTERVAL  OF  POSITION 

WRITEC  6 , 5 ) 

5  FORMATC/) 

WRITE( 6 , 31 )L 

31  FORMAT (5X. 'POSITION  STEP  NUMBER* ’,13) 

TSTEP=DX*DFLOAT ( L ) 

WRITEC  6 , 32 )TSTEP 

32  FORMATC 5X, ' X-POSITION* ' , FI  0 . 4/ ) 

WRITEC  6,101) 

101  FORMATC IX. ‘THE  U  VELOCITY  COMPONENT’) 

DO  300  KK=1 , NSTEP 
IS=CKK-1)*11+1 

IFIN=KK*10+KK 

WRITEC  6.82) CTHECI.L. 1 ) , I=IS  > IFIN) 

300  CONTINUE 

82  FORMATC 11C2X.F8. 5)) 

WRITEC  6,102) 

102  FORMATC IX, ’THE  V  VELOCITY  COMPONENT') 

DO  301  KK=1, NSTEP 

IS=CKK-1 )*11+1 
IFIN=KK*1 0  +KK 

WRITEC 6, 82 ) CTHECI.L, 2 ),I=IS, IFIN) 

301  CONTINUE 
WRITEC 6 ,103) 

103  FORMATC IX. 'THE  T  FIELD  VARIABLE  ') 

DO  302  KK=1, NSTEP 

IS=CKK- 1 )*1 1+1 
IFIN=KK*10+KK 


I 


,•»  ,♦»  »«  .<1 
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h 


i 


WRITE ( 6 ,82) (THE (I, L»  3) ,I=IS>IFIN) 
302  CONTINUE 
ITMAX=0 
20  CONTINUE 
RETURN 
END 
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